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INTRODUCTION 


Motivation  and  Overview 

The  difficulty  involved  in  experimentally  validating  NASP  type  high  altitude/high  Mach 
number  configurations  has  driven  a  large  hypersonic  CFD  effort.  The  nonequilibrium  real  gas 
effects  and  steep  gradients  found  in  this  regime  are  characterized  by  multiple  time  and  length 
scales.  Specifically,  the  relaxation  lengths  for  such  nonequilibrium  processes  often  differ  from 
the  convective  length  scale  by  several  orders  of  magnitude,  introducing  both  stiffness  and 
complexity  into  attempted  analysis  (1,31).  Moreover,  these  nonequilibrium  length  scales  may 
change  radically  throughout  the  field.  In  terms  of  numerical  simulation,  these  factors  drive  up 
the  cost  of  computation  by  grossly  increasing  the  resolution  requirements  for  acceptable 
solutions  (30,32). 

In  response  to  such  needs,  the  technique  of  adaptive  grid  embedding  locally  refines  the 
computational  mesh  based  on  information  from  a  developing  solution  (6,16,33).  Triggered 
by  both  gas  dynamic  and  nonequilibrium  variables,  embedding  increases  the  resolution  of 
flow  features  associated  with  these  processes.  In  this  way,  the  grid  scale  locally  adapts  to  the 
dominant  physical  scale  within  the  gas.  The  technique  ensures  grid  resolution  comparable 
with  all  important  physical  scales  of  the  problem  while  still  avoiding  unnecessary  resolution  in 
smoother  regions  of  the  flow. 

Of  course,  in  reacting  blunt  body  flows,  the  chemical  and  convective  lengths  are  not 
the  only  scales  in  the  flow.  The  bow  shock  adds  a  third  type  of  length  scale.  However,  since 
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this  structure  is  on  the  order  of  the  local  mean  free  path,  its  scale  is  too  disparate  to  attempt  to 
resolve.  In  fact,  with  the  exception  of  transatmospheric  problems,  the  bow  shock  may  be 
treated  as  a  discontinuity.  With  this  motivation,  shock  fitting  removes  the  shock  from  the 
domain  interior  in  the  present  work.  The  discontinuous  shock  assumption  is  consistent  with  a 
continuum  description  of  the  flow,  but  breaks  down  in  the  large  Knudsen  number  regime 
characteristic  of  very  high  altitude  flight  (28). 

From  a  computational  standpoint,  an  adaptive  embedding  procedure  requires 
unstructured  data  storage  and  a  compact  computational  stencil.  In  response  to  these 
pressures,  Ni's  (29)  node-based  finite-volume  scheme  was  chosen.  Over  the  course  of  the 
research  the  initial  in  viscid  perfect  gas  scheme  was  extended  to  include  chemical  source  terms, 
viscous  terms,  and  axisymmetric  configurations. 

Organization 

This  work  is  divided  into  three  main  parts.  The  first  takes  up  high  temperature  gas 
dynamics.  After  briefly  discussing  nonequilibrium  and  real  gas  phenomena  in  air,  it  develops 
the  governing  equations  and  gas  model  used  to  describe  chemically  reacting  high  temperature 
air.  The  model  applies  to  both  coupled  and  uncoupled  reacting  mixtures  and  degenerates  to  a 
classical  perfect  gas  when  the  chemistry  is  frozen.  The  final  results  of  these  discussions  are 
the  chemical  source  terms  for  use  in  the  governing  equations. 

In  Part  II,  focus  shifts  to  the  adaptive  numerical  algorithm  and  its  associated  numerical 
procedures.  The  integration  scheme  emphasizes  the  conservative  nature  of  the  governing 
equations.  Discussions  of  the  embedded  mesh  procedure  also  stress  this  point,  especially 
with  respect  to  the  treatment  of  computational  interfaces  between  embedded  mesh  levels.  This 


Part  closes  with  a  presentation  of  shock  fitting  algorithms  for  both  nonequilibrium  and 
equilibrium  real  gas  models. 

Part  III  contains  the  major  results  and  conclusions.  These  discussions  include 
analysis  of  both  fundamental  physical  phenomena  in  blunt  body  flows  and  an  overall  look  at 
the  effectiveness  of  the  adaptive  procedure  in  a  hypersonic  environment.  In  examining  the 
physical  process  found  in  high  temperature  shock  layers,  the  analysis  emphasizes  the  im¬ 
portance  of  length  scales  and  their  effect  on  the  flow's  behavior.  Such  investigations  gave 
experience  with  the  adaptive  method  as  an  engineering  tool,  and  they  provided  a  firm  basis  for 
evaluating  the  effectiveness  of  the  embedded  procedure. 
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I  REAL  GAS  DYNAMICS 


1 .  Some  Elements  of  High  Temperature  Gas  Dynamics 


We  briefly  discuss  some  aspects  of  high  temperature  gas  flows  and  present  an 
extended  form  of  the  conservation  equations  to  model  such  flows.  This  permits  close  study 
of  the  terms  responsible  for  modeling  the  physics  of  gases  at  elevated  temperatures.  The 
examination  reveals  the  strengths  and  weaknesses  in  the  present  modeling,  and  provides 
insight  into  which  physical  processes  may  be  described  within  this  framework. 

Ehysical  Processes 

Nonequilibrium  involves  both  the  chemical  composition  and  the  internal  energy 
storage  modes  of  the  gas  particles.  In  steady,  air  breathing  flight  within  the  Earth's 
atmosphere,  hypersonic  vehicles  experience  both  kinds  of  nonequilibrium.  Figure  1.1  shows 
the  corridor  for  sustained  flight  with  regions  of  thermo-chemical  excitation  superimposed 
(20).  Despite  the  somewhat  approximate  boundaries  in  the  sketch,  it  is  clear  that  these  modes 
affect  nearly  all  flight  at  hypersonic  Mach  numbers. 

Nonequilibrium  Rate  Processes 

While  excitation  of  internal  modes  is  important,  excitation  alone  does  not  imply  a 
relevant  nonequilibrium  process.  In  the  context  of  flowing  gases,  nonequilibrium  implies  that 
the  time  for  an  energy  mode  to  equilibrate  is  an  appreciable  fraction  of  some  time  scale  within 
the  flow.  On  average,  chemical  reactions  require  thousands  of  collisions  between  the  re¬ 
actants.  These  collisions  occur  over  a  finite  time  and  distance.  In  general,  the  more  efficient 
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energy  transfer  between  internal  modes  of  a  particle  implies  less  space  and  time  to  achieve 
equilibrium.  This  observation  often  leads  to  a  heavy  particle  temperature  approximation  (21). 
The  temperature  characterizes  a  particle's  internal  energy  by  assuming  thermal  equilibrium 
takes  place  instantaneously.  That  is,  all  internal  modes  remain  in  equilibrium  with  each  other. 
While  valid  throughout  much  of  the  flight  corridor  in  Figure  1.1,  recent  work  by  Candler  (5), 
Park  (31),  and  others  (4)  suggests  that  the  approximation  breaks  down  near  orbital 
conditions. 


Mach  Number 

FIGURE  1.1 

Velocity  -  Altitude  map  with  legions  of  modal  excitation  superimposed  on  it.  (From  Ref.  20.) 


At  elevated  Mach  numbers  and  low  densities,  the  vibrational  temperature  may  take  appreciable 
time  (and  distance)  to  equilibrate  with  the  translational-rotational  temperature.  This  results  in 
thermal  nonequilibrium  within  the  different  molecular  species  in  heated  regions  of  the  shock 
layer.  In  this  context,  "molecular  species"  refers  mainly  to  nitrogen,  since  oxygen  molecules 
tend  to  dissociate  before  their  vibrational  and  translational  temperatures  differ  appreciably.  It 
is  worth  noting,  however,  that  the  stagnation  enthalpy  error  produced  by  ignoring  this  tem¬ 
perature  difference  is  negligibly  small  for  a  wide  class  of  problems.  Thus,  while  multiple 
vibrational  temperatures  remain  important  in  re-entry  or  transatmospheric  problems,  a 
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one-temperature  model  does  accurately  model  a  wide  class  of  problems  within  the  Earth's 
atmosphere. 


1.1  Conservation  Equations  for  Chemically  Reacting  Viscous  Flow 

Many  texts  develop  the  conservation  equations  for  chemically  reacting  viscous  flows 
(for  example  (36),  (13).  (1)).  This  section  first  states  this  set  and  then  discusses  specific 
modeling  for  terms  within  these  relations. 


Global  conservation  of  mass  and  momentum  simply  account  for  mass  and  momentum 
flux  through  a  fluid  without  concern  for  any  internal  structure  of  the  gas  particles.  The 
description  therefore  remains  unchanged  from  the  classical  Navier-Stokes  or  Euler 
formulation.  In  full  conservation  form: 


Continuity 

^  +  V.(pV)  =  0 


(1.1) 


x  Momentum 


d(pu)  3 jpu2  +  p)  djpuv)  _  B(rxx)  d(TXy) 

df  dx  By  dx  By 


y  Momentum 


3(pv)  Bjpuv)  B{pyi  +  p)_B( ryx)  B{ ryy) 
Bt  Bx  By  Bx  By 


(1.2) 


Here,  the  shear  stresses  are 


and  the  normal  viscous  stresses  are 
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The  Stokes  hypothesis  eliminates  the  second  coefficient  of  viscosity  by  setting  A  = 


iyy 


-141-1) 


Since  the  energy  equation  imposes  a  detailed  balance  of  all  energy  fluxes,  it  must  include  an 
explicit  description  of  the  diffusive  energy  flux  arising  from  gradients  in  species 
concentrations.  However,  collecting  those  effects  within  the  heat  flux  term,  q,  the  equation  to 
be  couched  in  classical  form. 


f +a=£d+  *£*-5^  +  '*-'>*&*  +  -  *> 


(1.3) 


Here  e  refers  to  the  total  internal  energy  per  unit  mass,  and  the  stagnation  enthalpy 
becomes: 

P  (1.4) 


As  mentioned,  q  includes  both  the  Fourier  heat  conduction  term,  and  the  diffusive 
energy  flux  contributions. 

q  =  -kVT+J\  piVi hi 

s  (1.5) 


Here,  U,-  is  the  diffusion  velocity  of  the  <th  component  in  the  mixture.  More 
precisely,  U measures  the  relative  velocity  of  the  i*  species  to  the  bulk  velocity  of  the  gas. 

The  diffusion  velocity  also  contributes  to  the  species  conservation  equations.  These 
relations  trace  the  net  production  or  destruction  of  each  species  in  the  mixture. 
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For  the  1th  species: 

d(P«)  ,  %P«“)  3(p,v)  _  w,  3(p,g,)  3(p«v.) 

9r  9*  "  4  dx  dy  (L6) 

The  source  term  Wi  contains  the  net  reaction  rate  expression,  and  in  general  is  a  summation 
over  all  reactions  contributing  to  the  net  rate  of  change  of  any  particular  species. 


Naturally,  conservation  of  atomic  nuclei  and  electrons  prohibits  the  gross  numbers  of 
particles  of  any  species  from  changing. 


=0 


(1.7) 


Additionally,  global  continuity  forces  the  net  diffusive  flux  of  all  individual  species  to 


cancel. 


X  pc, Ui  =  0 


(1.8) 


State  Vector  Form 

The  conservative  form  of  the  governing  equations  permits  an  arrangement  more 
amenable  to  numerical  computation. 

_9R  as  w 

dr  dx  dy  "  dx  dy  dm 


Here 


U  = 


pu 

pv 

p 

pu 2  +  p 

puv 

pu 

pv 

puv 

pv 2  +  p 

e 

,  F  = 

u(e+p) 

,  G  = 

v(e+p) 

pi 

piU 

piv 

_  Ps  _ 

• 

• 

• 

_  P*“  _ 

_  P'v  _ 
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and 


0 

r  0 

'  0  ' 

Ixx 

0 

*xy 

.  T» 

0 

R  = 

(urn  +  Vfy  -  qx) 

,  S  = 

(u*xy  +  VTyy-  qy) 

,  w  = 

0 

piUi 

piji 

Wi 

Ps% 

PsVs 

.  Ws  . 

F  and  G  are  the  inviscid  flux  vectors  while  R  and  S  contain  the  viscous  fluxes.  Setting  the 
viscous  flux  vectors  identically  to  zero  results  in  the  governing  equations  for  chemically 
reacting  inviscid  gas  mixtures.  W  contains  the  chemical  source  terms  for  either  case.  Setting 
this  vector  to  zero  degenerates  the  set  of  equations  further  and  results  in  the  Euler  equations 
for  perfect  gas  modeling. 

Normalization 

Table  1.1  shows  the  scaling  of  flow  variables  and  properties  used  to  non-dimension- 
alize  the  governing  equations.  The  subscript  (  )f  refers  to  frozen  conditions. 


TABLE  1.1 

Scaling  of  flow  variables  and  properties  used  to  non-dimensionalize  the  governing  equations. 


Reference 

quantity 

HHBKSSSZSSI&ISHHIH 

p 

P- 

i 

U,  V 

°Foo 

A/ x<x>  M y— 

e 

Dependent  upon  modeling 

P 

T 

1 

t 

/afot 

— 

P 

P~ 

1 

k 

L. 

1 

Dim 

Dim* 

1 

X 

Rn 

_ 

y 

Rn 

— 

Taking  the  frozen  sound  speed  as  a  reference  velocity  conveniently  rescales  the  freestream 
flow  velocity  and  Mach  number  to  numerically  equal  values. 


Non-Piroensionalization  of  Diffusion  Properties 

Inserting  this  scaling  into  the  governing  equations  permits  a  collection  of  all 
dimensionless  parameters  within  R,  S,  and  W.  Collecting  these  parameters  into  the  non- 
dimensional  heat  flux  vector  q  and  shear  stress  tensor  f  leaves  the  form  of  the  governing 
equation  set  (1.9)  unchanged. 


Defining  the  Reynolds,  Prandtl,  and  Schmidt  number  (respectively)  as 

/X.  k. 


Sc  s- 


P- 


P-D i 


unm 


permits  r  to  be  expressed  without  dimension. 


(1.10) 


xyy  ~  “  u*) 

Txy  =  +  v;) 


where,  H*  =  u*  =  and  v*  =  . 


Since  significant  ionization  occurs  only  at  temperatures  too  high  to  simulate  accurately 
within  the  present  single  temperature  framework,  any  charge  separation  predicted  due  to  ionic 
species  must  remain  small.  Given  this,  Fick's  first  law  of  diffusion  relates  the  diffusion 


velocity  to  species  gradients. 


C,  = 

Pi 


(1.11) 


Collecting  dimensional  parameters  in  the  energy  equation  returns  a  non-dimensional  form 
of  this  quantity.  In  cartesian  form  the  diffusion  velocities  become: 


‘  teSc\  mp*dxmJ  '  ReSc[  'mp'idy' 


(1.12) 


With  these  quantities  now  known,  we  may  write  the  non-dimensional  heat  flux  vector. 


q\PrRe^WT 


Breaking  this  vector  into  its  cartesian  components,  q  becomes; 
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(1.13) 


(  3T  _  v1  phjfiimdCjX 

Qx=  \Pr  Re  yp~  dx  T  *<?  Sc  dx  ) 

(  dT  y  phiDimdCj\ 

qy~\PrRe^dy  Re  Sc  By J 


After  developing  a  mixture  gas  model,  the  next  chapter  details  the  non-dimension- 
alizadon  of  chemical  source  terms. 


Precise  modeling  of  diffusion  properties  requires  details  of  the  collision  cross-sections 
from  kinetic  theory  (21,  36).  However,  the  level  of  accuracy  for  the  current  modeling  does 
not  warrant  such  details.  Moreover,  fj.  and  k  should  take  into  account  mixture  values  in 
reacting  calculations,  but  such  effects  are  of  higher  order,  and  are  not  crucial  for 
understanding  the  basic  physics  here. 


Sutherland's  law  predicts  fi  to  within  10%  below  9000  K  (10).  All  viscous 
calculations  presented  later  use  this  approximation. 


*_  P  JT-  +  HO.41/ T  \| 
1  “/u=[r+  no.4  J\rJ 


(1.14) 


The  definition  of  Prandtl  number  relates  thermal  conductivity  to  viscosity. 

^  _  P~cprm 
Pr m 


A  constant  Prandtl  number  assumption  provides  a  first  approximation  for  the  high 
temperature  transport  property  behavior, 


!c*  53  k  -  CP'-  ~  „*  cp 
k-  hL  cp- 
Prm 


(1.15) 
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Unless  otherwise  stated,  all  viscous  calculations  assume  that  Cp/cPFj=  1.0. 


Similarly,  the  definition  of  Sc  provides  a  basis  for  evaluating  Dim 


t 


p~Sc 


' im 


for  constant  Schmidt  number  modeling. 


(1.16) 


The  Degree  of  Noneouilibrium 

The  earlier  statement  that  nonequilibrium  is  relevant  when  relaxation  occurs  over  an 
"appreciable"  distance  was  accurate,  but  not  very  quantitative.  A  clearer  statement  of  the 
degree  of  nonequilibrium  relates  the  fluid  and  chemical  time  scales 


Consider  a  steady  state  inviscid  species  equation.  From  Equation  (1.6), 


(1.17) 


This  may  be  non-dimensionalized  as  before  with  freestream  conditions. 

ai£iiir ,  RjrtVj a.i8) 

dX 

Here,  IP  is  a  form  of  the  Damkdhler  number  and  relates  the  chemical  reaction  rate  to 
the  bulk  fluid  motion.  The  reference  reaction  rate,  Wjjef,  was  evaluated  at  conditions 
downstream  of  a  normal  shock.  As  *P  increases,  the  chemical  length  scale  shrinks  until 
identical  equilibrium  modeling  describes  the  flow.  Similarly  as  *P  approaches  zero  the  flow 
becomes  progressively  more  frozen. 
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1.3  Chemical  Effects  in  High  Temperature  Air 


With  the  governing  equations  for  a  general  nonequilibrium  chemically  reacting  gas 
mixture  completely  outlined,  we  next  exanine  how  they  model  a  known  system.  Moreover, 
the  chemical  composition  of  high  temperature  air  will  provide  a  datum  by  which  to  judge  the 
results  of  later  computations.  For  equilibrium  air,  the  species  equations  reduce  to  the 
equivalent  Law  of  Mass  Action  form. 


Temperature  K 

FIGURE  1.2 

Compcwition  of  equilibrium  air  u  a  function  of  temperature  at  1  atm  pressure  (1). 


Figure  1.2  illustrates  the  temperature  dependence  of  the  equilibrium  mole  fraction  X  for  the 
five  major  constituents  of  air.  The  mole  fraction  of  the  i*  species  1 f/  is  defined  as  the  ratio  of 
the  number  of  moles  of  i,  to  the  total  number  of  moles  in  the  mixture  JS. 

At  atmospheric  pressure,  the  system  remains  quiescent  from  a  few  degrees  above 
absolute  zero  to  nearly  800  K.  Over  this  range,  rotation  and  translation  are  fully  excited,  and 
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the  equipartition  principle  dictates  that  both  modes  participate  in  energy  storage;  At  about 
800  K,  the  vibrational  mode  begins  to  contribute  to  such  storage,  but  overall  energy  levels  are 
insufficient  to  appeciably  dissociate  any  of  the  vibrating  molecules.  Classically  speaking,  this 
vibrational  excitation  alters  the  specific  heat  ratio  •  affecting  such  things  as  the  relationship 
between  pressure  and  internal  energy  of  a  static  gas. 

At  about  2000  K,  these  relatively  benign  effects  change  character.  Molecular  vibration 
and  collisions  have  now  increased  to  the  point  where  appreciable  oxygen  dissociates.  Again 
the  specific  heat  ratio  begins  to  change  rapidly,  except  that  now  it  varies  more  wildly,  since 
changes  in  temperature  also  bring  about  changes  in  composition. 

By  4000  K,  only  trace  amounts  of  O2  remain,  and  N2  begins  to  dissociate.  This 
process  continues  until  essentially  only  atomic  species  survive  above  9000  K.  At  such 
temperatures,  the  equilibrium  constants  favor  atoms,  and  the  diatomic  species  which  do  form 
exist  only  long  enough  to  dissociate  again.  Above  this  temperature,  additional  thermal  energy 
can  only  be  stored  by  breaking  down  the  electronic  bonds  holding  electrons  to  the  atomic 
nuclei. 


One  interesting  feature  is  the  bubble  of  NO  appearing  between  2000  and  6000  K.  As 
background,  consider  the  five  basic  neutral  species  reactions  in  the  air  system  (below 
9000  K). 


i. 

n2  +  m 

2N  +  M 

2. 

02  +  M 

20 +  M 

3. 

NO  +  M 

«=*  N  +  O  +  M 

4. 

NO  +  O 

c=> 

N  +  02 

5. 

0  +  N2 

<=> 

NO  +  N 
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The  last  two  exchange  reactions  (or  "shuffle"  reactions)  are  responsible  for  creating 
NO  in  the  system.  Atomic  oxygen  collides  with  nitrogen  molecules,  forming  NO  which 
dissociates  more  slowly  than  it  is  produced,  shifting  the  equilibrium  composition  upward. 
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FIGURE  t-J 

Equilibrium  constants  for  neutral  species  of  air.  (From  Ref.  33.) 

Figure  1.3  provides  greater  insight  into  the  role  of  the  exchange  reactions.  This  shows  the 
equilibrium  constants  resulting  from  the  ratio  of  the  net  forward  to  the  backward  rates  of  the 
five  reaction  system.  Note  that  at  relatively  lower  temperatures,  the  shuffle  reactions  proceed 
at  rates  orders  of  magnitude  faster  than  the  dissociation  paths.  At  those  temperatures,  the 
reverse  rates  tend  to  dominate  both  (4)  and  (5),  and  the  NO  produced  by  (4)  reacts  in  (5), 
tending  to  hold  the  overall  level  of  NO  nearly  constant.  Additionally,  both  reactions  consume 
atomic  nitrogen,  further  depleting  whatever  trace  amounts  may  exist. 

Although  completed  at  standard  pressure,  the  composition  plot  (Fig  1.2)  holds 
qualitatively  to  quite  low  pressures.  Lowering  the  pressure  shifts  the  curves  to  the  left, 
raising  the  degree  of  dissociation  at  any  fixed  temperature.  Finally,  this  figure  provides  some 
reference  state  for  the  high  temperature  shock  layers  presented  later. 
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2.  Modeling  High  Temperature  Air 


The  gas  model  developed  to  describe  high  temperature  reacting  flow  extends  the  concept 
of  a  finite  rate  ideal  dissociating  gas  (1 1, 23)  to  include  both  coupled  and  uncoupled  systems. 
While  capable  of  predicting  the  essential  physical  phenomena,  this  model  introduces  a 
minimum  of  algebraic  complexity.  Strictly  speaking,  it  is  reasonably  accurate  to  9000  K,  and 
presently  includes  only  electrically  neutral  species.  Nevertheless,  Chapter  6  presents  evidence 
which  demonstrates  accurate  modeling  of  flows  where  small  portions  of  the  shock  layer 
include  much  higher  temperatures. 

After  examining  mixture  thermodynamics,  attention  turns  to  the  formulation  of  chemical 
source  terms  and  details  of  the  chemical  modeling.  Although  the  focus  is  on  air  chemistry,  an 
extension  to  general  mixtures  is  straightforward. 


2.1  Thermodynamics  of  a  Mixture,  of  Ideal  Dissociating  Gases 

The  relationship  between  thermodynamic  quantities  in  a  system  relies  on  the  details  of  gas 
behavior  on  a  molecular  scale.  The  primary  simplifications  in  the  present  modeling  concern 
internal  energy  storage  within  molecular  species.  A  single  temperature  characterizes  all  internal 
modes  since  the  vibrational  mode  follows  Lighthill's  "half  excited"  assumption.  In  1983,  Park 
et  al  demonstrated  that  this  assumption  accurately  predicts  static  enthalpy  over  the  range  of  tenv 
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peratures  where  vibrational  excitation  is  present  (30).  Figure  2.1  reproduces  this  work, 
comparing  the  static  enthalpy  and  vibrational  energy  of  Lighthill's  model  with  a  more  accurate 
modeling  (14).  The  indication  is  that  only  a  small  error  exists  over  much  of  the  shock  layer. 


FIGURE  2.1 

Energy  level  companion  of  equilibrium  air  between  half-excited  model  and  more  precise 
model  of  reference  (14)  (reproduced  from  ref.  (31)). 


Pressure 

For  ideal  gas  mixtures,  Dalton's  Law  of  Partial  Pressures  states  that  the  total  pressure  is 
the  sum  of  the  contributions  of  each  component 
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Here,  the  density  of  the  /th  species,  p, ,  is  related  to  the  total  density  by 


(2.2) 


Since  the  law  is  purely  mechanical,  it  remains  unchanged  in  nonequilibrium  flows. 

Internal  Energy  and  Enthalpy 

The  overall  internal  energy ,£,  for  an  inert  mixture  is  the  sum  of  the  component  internal 
energies. 

£  »  X  £, 

s 

Here  the  tilde  indicates  that  reactions  are  absent.  The  mixture  enthalpy  is 


making  the  total  internal  energy  per  unit  volume  of  a  mixture 


e 


With  this  definition,  the  total  enthalpy  becomes 


(2.3) 


(2.4) 


If,  however,  the  mixture  is  chemically  reacting,  the  internal  energy,  E,  must  track  the  energy 
absorbed  during  dissociation.  That  is: 


In  which  &dj  is  the  characteristic  temperature  for  dissociation  of  they*  molecular  species. 

Equation  (2.5)  implies  that  a  mixture  of  only  atoms  will  have  zero  energy  at  zero  Kelvin. 
While  convenient  for  some  applications,  this  datum  is  not  unique.  Moreover,  since  at  low 
temperatures  air  is  virtually  all  molecules,  we  reset  the  datum  to  zero  internal  energy  at  absolute 
zero  and  100%  molecules. 

s  J-  "$s£r  (2.6) 

Energy  Storage 

Expressions  for  the  internal  energy  of  each  species,  £/,  stem  from  the  kinetic  theory  of 
gases.  Each  species  stores  ^RfT  for  each  equilibrium  degree  of  freedom. 

£ Trans  =  2  »  Eftot  ~  ,  and  Evib  =  ra^E 

However,  since  the  current  model  assumes  vibration  is  only  half  excited,  Evib=  \rA2t- 

In  a  system  with  atomic  species.  A,  and  molecular  species,  A2 

3 

Atomic  Species:  E  =  Efrons  ~  2  ^A^ 

3  1 

Molecular  Species:£  =  Efrons  +  Er01  +  Evib  =  5  R*2T  +  ^2^  +  2  ^A2T  =  ^  Ra2T 


For  the  specific  mixture  of  the  five  major  species  in  high  temperature  air  N,  O,  NO,  N2,  and 
O2  (species  1-5  respectively),  the  internal  energy  becomes 


As  a  check,  consider  the  special  case  of  a  symmetric  dissociating  gas,  A2  <=>  2A,  with  charac¬ 
teristic  temperature  &d  (i«-  with  C2  =  03  =  C5  =  0,  C4  =  1  -  c\,  and  m 4  =  2mi), 

£=/?„, (37- +c„eD) 
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This  is  precisely  Lighthill’s  result  for  the  same  case  (36). 


Finally,  the  integration  scheme  uses  the  total  internal  energy  per  unit  volume,  based  on 


combining  (2.3)  and  (2.7) 


(2.8) 


2.2  Law  of  Mass  Action 

The  "half  excited"  assumption  corresponds  to  certain  behavior  of  the  molecular  partition 
function.  This  is  expected,  since  those  functions  statistically  describe  how  the  energy  storage 
modes  of  a  particle  behave.  However,  Lighthill's  original  work  considered  only  symmetric 
diatomic  gases,  and  a  model  of  the  complete  air  system  must  extend  the  analysis  to  include 
asymmetric  particles  such  as  NO. 

For  the  general  dissociation  reaction,  AB  <=>  A  +  B,  the  number  densities  measure 

Nab  =  Total  number  of  nuclei  per  unit  volume  (including  those  in  both 
atoms  and  molecules) 

NA  s  Number  of  A  atoms  per  unit  volume 

NB  s  Number  of  B  atoms  per  unit  volume 

NAB  s  Number  of  AB  molecules  per  unit  volume 

Within  this  present  framework,  the  total  number  of  nuclei  remains  constant 

2 Nab  +  Na  +  Nb  =  Nab  (2.9) 


21 


From  statistical  mechanics,  the  law  of  mass  action  for  this  general  system  relates  the  number 
densities  to  the  partition  functions  of  each  species  and  the  change  in  zero  point  energy,  of 
the  reaction. 


NANB  ^Q*Qb 
Nab  Q^b 


(2.10) 


The  change  in  zero  point  energy  divided  by  Boltzman's  constant,  K,  is  the  characteristic  tem¬ 
perature  for  dissociation,  referred  to  earlier  as  6fo. 

We  may  now  define  the  degree  of  dissociation  as  the  ratio  of  particles  of  a  particular  type  to  the 
total  number  of  nuclei  available. 


aAs 


(2.11) 


When  combined  with  a  statement  of  nuclei  conservation, 

NAB  *  (1  -aA~  ocb)~^- 

we  may  define  O-ab  as 


(2.12) 


(2.13) 


Note  that  in  single  reaction  systems  with  stoichiometric  coefficients,  little  distinction  is  usually 
made  between  the  degree  of  dissociation,  a,-,  and  the  mass  fractions,  c,\  However,  if  other 
kinds  of  particles  are  present,  the  degree  of  dissociation  differs  from  the  mass  fraction,  as  will 
be  detailed.  For  now,  simply  note  that  by  definition  aA  +  olb  +  <*AB  -  but  cA  +  cb  +  cab 
is  not  unity  in  the  presence  of  other  particles. 

Repeating  the  law  of  mass  action  in  terms  of  the  degree  of  dissociation  and  the  total 
number  of  nuclei.  Nab  . 
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*kCtB  _  Q^Q8  C-Ab It 
<*ab  2  NABQ*b 


(2.14) 


These  nuclei  have  an  average  mass  defined  by 


,  _  (mA  +  ms)NAB  +mANA  +  mBNB 
Nab 


(2.15) 


m'  remains  constant  for  any  reaction  at  any  specific  time  provided  that  the  reactants  exist  at 
stoichiometric  ratios.  Moreover,  for  any  reaction  in  a  multi-reaction  system,  we  may  define  a 
density,  p\  derived  from  the  fraction  of  the  mixture's  density  that  is  actually  involved  in  any 
specific  reaction. 


P' s  P(C/t  +  cb  +  cab) 


(2.16) 


From  these  above  two  definitions, 


(2.17) 


and  the  law  of  mass  action  becomes 


-  m'  Q*QBe-aJT 
<Xab  2 p'V  qab 


(2.18) 


Notice  again  that  since  (Xa,  as,  and  (Xab  represent  the  degree  of  dissociation,  this  form  applies 
for  systems  containing  other  (possibly  inert)  reactions. 


In  the  very  special  case  of  asymmetric  dissociating  reaction  with  no  other  species  present 
and  reactants  at  stoichiometric  ratios,  Equation  (2.18)  becomes 


a 2 

1  -  a  \2V  qAA  I  p 


(2.19) 


The  contents  of  the  parenthesis  on  the  right  side  have  units  of  mass  per  volume.  Moreover,  for 
nitrogen  and  oxygen  systems,  the  ratio  of  partition  functions  remains  remarkably  constant 
below  9000  K.  By  defining  a  characteristic  density ,  po ,  we  may  approximate  this  term  with  a 
constant,  thus  avoiding  the  added  complexity  of  evaluating  partition  functions.  It  is  easy  to 
show  (36,  pp.  159)  that  this  simplification  is  equivalent  to  the  "half  excited"  assumption  stated 
earlier. 


W  vPjfiTTrauWi  W  t*5n  Hi 


Returning  to  the  boxed  Equation  (2.18),  we  seek  to  define  a  characteristic  density  for  this 
asymmetric  reaction.  However,  the  form  of  (2.18)  is  unfortunate  since  m'  will  vary  when 
other  reactions  are  present.  Specifically,  in  a  general  mixture  nothing  insures  that  a*  will 
equal  ,  and  thus  there  is  the  possibility  of  an  excess  of  particles  which  will  remain  inert. 
Any  attempt  to  define  a  po  parameter  will  result  in  a  quantity  which  varies  according  to  m. 

In  order  to  reexpress  —  examine  the  definition  of  p'  (2. 16).  Recall  the  mass  fraction  of 
the  1th  species  is  the  total  mass  of  i  divided  by  the  total  mass  in  the  system. 


P’  *  p(ca  +  cB  +  cAb),  but  c,  s 


mjN 1 

X  miNi 


where 


Substituting  for  the  mass  fractions, 


P’  =  P 


thaN*  +  msNB  +  mAaNAB 

X  mjN* 


(2.20) 


in  a  mixture  with  s  components.  The  total  number  of  nuclei  in  this  system  stems  from 
Equation  (2.9): 


24 


Na.  -  I  m,N'  +  m,N‘  +  %%  m,N‘ 


(2.21) 


We  now  re-form  the  ratio  ~  from  (2.10)  and  (2.16). 


(2.22) 


P’  + 

F\mA  mB  mAB 


Defining  m  as  the  average  mass  per  nucleus  in  AB  molecules, 

n  =  mAt.me  SM 


and  multiplying  both  the  numerator  and  denominator  of  (2.22)  by  this  quantity  results  in 


ml  =  mllm)  _ _ m _ _ _ =  m. 

p'  ^ m  ? 


(2.23) 


m  is  constant  for  any  reaction  irrespective  of  any  other  inert  or  active  species  present.  Thus, 
the  variation  has  moved  from  m'  into  ~p,  while  simultaneously  re-expressing  the  ratio  of 
(2.23)  in  terms  of  the  (more  common)  mass  fractions. 

To  demonstrate  the  value  of  this  step,  re-examine  the  law  of  mass  action.  Equation 


(2.18)  becomes: 


aA^L  - 1  m  Q^QB  e-eJT 
ccab  p\2V  qAB  J 


(2.24) 


Following  Lighthill,  the  term  in  the  bracket  may  be  defined  as  a  characteristic  density  for  disso¬ 
ciation  of  AB,  poAB- 


ym&Q8 

2V  qAB 


(2.25) 


This  form  reduces  identically  to  that  presented  by  Lighthill  in  1957  for  a  symmetric  diatomic 
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Values  of  On  and  Bn 

Figure  2.2  displays  plots  of  Equation  (2.25)  for  the  major  diatomic  species  in  equilibrium 
air  below  9000  K.  Reference  36  provided  the  values  used  in  evaluation  of  the  partition 
functions. 


♦*  Niiro*en  Oxide  0-  Oxygen  ■»  Nitrogen 


Temperature  K 


FIGURE  2.2 


Characteristic  density .pQ  for  the  three  major  diatomic  molecules  in  high  temperature  air. 


Table  2. 1  summarizes  values  of  po  ,  and  Bd  for  N2, 0 2,  and  NO. 

TABLE  2.1 

Characteristic  temperatures  and  densities  for  dissociation  of  the  major  molecular  species  in  equi¬ 
librium  air  below  9000  K. 


j  A/2 

- & - 

NO 

Characteristic  Temp.  (K) 
Characteristic  Density  (kg/m^) 

113000 

130000 

59  500 

150  000 

75  500 

30  000 
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2.3  Finite  Rate  Expressions 


Symmetric  Diatomic  Gas 

Freeman  (11)  approximated  the  general  chemical  rate  expression  by  separating  the  law  of 
mass  action  into  forward  and  backward  components  (37,  pp.  232).  For  an  ideal  dissociating 
gas  with  law  of  mass  action 

a 2  -  P°  e-edr 

l~a  P  (2.26) 


the  rate  of  dissociation  da/d,  is: 


=  CTiplll  -  a)e-^T  _  P-a 2 
dr  r '  on 


(2.27) 


which  describes  the  production  of  atoms.  Notice  that  for  symmetric  diatomic  gases  with  no 
other  species  in  the  mixture,  a  -  ca- 


Mixture  of  Ideal  Dissociating  Cases 

The  law  of  mass  action  for  asymmetric  molecular  dissociation  with  a  half  excited 
vibrational  state  (2.24)  may  be  re-written  in  terms  of  the  characteristic  density  and  degree  of 
dissociation. 

=gg.e-ft/r 

ttAB  P  (2.28) 


The  analog  of  (2.27)  is  simply 


iaAB 

dr 


=  c™ 


-2-aAaB  -  aABe 
Pd 


(2.29) 


which  describes  the  production  of  AB  molecules.  Note  that  p  =  Zj  Pi  as  in  (2.2). 
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Equation  (2.29)  describes  AB  production  as  a  function  of  the  degree  of  dissociation 
instead  of  mass  fraction.  Since  the  species  equations  express  W  in  terms  of  the  mass  fraction, 
(2.29)  requires  reformulation. 

Recalling  the  definition  of  ocab  from  (2.13)  and  then  substituting  for  N‘  and  Nab  leads  to 

daAB  _ _ 1 _ d Cab 

to  [C*B+ 9]  to 

Re-arranging  to  express  the  mass  fraction  in  terms  of  the  degree  of  dissociation  gives 

d cab  -PAQab 

to  P  to  (2.30) 

after  substituting  for  ~p  from  (2.23). 

For  the  symmetric  diatomic  case  ^  reduces  to  unity,  ca  +  cb  +  cab  -  1,  and  mA  -  mB. 
Equation  (2.30)  collapses  to  the  rate  expression  (2.27). 


Finally,  for  the  general  dissociative  reaction,  AB  <=>  A  +  B,  in  a  mixture  with  many 
components 

^AB-^CT^p 


^-aAaB  -  ciABe'^h 
Pd 


(2.3H 


This  simple  form  appears  to  rely  on  fortuitous  cancellation  and  careful  definition  of  p  .  How¬ 
ever,  by  expressing  the  degree  of  dissociation  in  terms  of  mass  fractions  through  the  number 
densities  M  and  total  number  of  nuclei  Nab>  one  may  show 

=  and  = 

P  AP  m“P  (2.32) 


Substituting  back  into  (2.31)  expresses  the  rate  of  change  of  species  mass  fractions  as  direct 
functions  of  these  mass  fractions  of  the  reactants  and  products. 
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(2.33) 


^£M.~-CTnp 


cab^ 


mAmBlpD 


CaCb 


Here,  again,  m  is  the  average  mass  per  atom  in  the  AB  molecule. 

Elimination  of  Species  Equations 

With  the  species  production  terms  formulated,  attention  returns  to  the  governing 
equations.  In  general,  as  many  species  equations  exist  as  species.  However,  for  a  system 
with  only  two  types  of  nuclei,  conservation  permits  discarding  two  differential  equations  - 
slightly  reducing  computational  requirements  and  effort.  In  the  present  work,  with  the  species 
N,  O,  NO,  N2  and  O2  numbered  1-5  (respectively),  we  elect  to  eliminate  the  last  two  species 
based  on  conservation  of  atoms. 


The  simplification  is  a  convenience.  Moreover,  since  we  presently  consider  only  the  air 
system,  the  notation  will  refer  to  N  and  O  directly,  dropping  the  pretense  of  A  and  B. 


Nitrogen  nuclei:  Nn  =  2 NNl  +  Nn  +  Nno 
Oxygen  nuclei:  No  =  2 N°*  +  N°  +  Nno 
Total  number  of  nuclei:  Nno  =  Nn  + No  (2.34) 

Differentiating  and  re-organizing  these  equations  leads  directly  to  the  rate  expressions  for 
molecular  nitrogen  and  oxygen. 

dcy2  _  (dcy  +  my  dfivoj 

dr  'dr  mno  dr  »  (2  35) 

d£Oz  _  (dcp  mp  dc^p] 

dr  \  dr  mNO  dr  ' 


In  addition  to  the  rate  expressions,  the  concentrations  of  N2  and  O2  need  to  be  expressed  as 
functions  of  the  other  species.  Taking  air  as  79%  nitrogen  and  21%  oxygen  by  number  (i.e. 
Nn  =  0.79 Nno,  Nq=  0.21Nno )  conservation  of  atomic  nuclei  yields 


cn2 


' _ 0-79 

0.79  mN  +  Q.2\mo 


cno 

mNo 


mN-CN 


(2.36) 
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0l  Io.79m*  +  O.21m0  "Wof^  0 


Non-dimensionalization  of  the  source  expressions  follows  directly  from  the  assumed 
scaling  (Table  1.1).  In  a  mixture  with  only  one  reaction,  the  inviscid  species  conservation 
equation  in  one  dimension  has  the  form 


^-  +  u~-  =  ^rv{[  F] -[  B  ]}  = 

3 1  Bx  mi  1  1  p 


(2.37) 


in  which  [F]  and  [B]  are  forward  and  backward  (i.e.:  dissociation  and  recombinat  ■>u/ 
dimensionless  expression  terms,  as  in  Eq.  (2.33),  and  C  =  C//nti  as  shown  in  Ref.  (36). 

After  rescaling  this  equation  with  the  parameters  in  Table  1.1,  the  source  term  becomes: 


K,.=  <frV([F]-[B])r 
'  V 


(2.38) 


Here,  the  star  denotes  non-dimensional  quantities  and  <t>  is  a  non-dimensional  rate  parameter 


defined  by : 


<j>(.  =  Cf‘  TrefPr‘f  W 

miUref 


(2.39) 


The  appearance  of  m  in  this  expression  is  consistent  with  the  reference  mass  appearing  in  the 
forward  and  reverse  terms.  The  generalization  extends  the  classical  definition  of  0  to  include 
systems  with  multiple,  coupled  reactions. 
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2.4  Chemical  Source  Terms 


Reactions  Considered 

The  chemical  source  terms  for  any  particular  species  results  from  the  summation  of  all 
rate  expressions  contributing  to  that  species.  The  primary  neutral  reactions  occurring  in  air 
below  9,000  K  are: 


1) 

n2  +  m 

0 

2N  +  M 

2) 

02  +  M 

<=> 

20 +  M 

3) 

NO  +  M 

<=> 

N+O  +  M 

4) 

N2  +  0 

<=> 

NO  +  N 

5) 

NO  +  0 

<=> 

02  +  N 

This  system  neglects  ions  of  NO  which  may  appear  at  such  temperatures  in  small  quantities. 
These  ions,  and  their  associated  free  electrons,  however,  contribute  primarily  to  charge 
separation  and  electron  densities  within  the  shock  layer,  and  are  assumed  to  be  of  secondary 
importance  in  the  present  study,  which  therefore  includes  only  the  neutral  reactions  listed 
above. 


Complete  RateJasBESSsiaps 

The  rate  expression  for  each  reaction  makes  use  of  the  species  production  terms  of  the 
previous  section  (e.g.  Eq.  2.35).  The  numerical  subscripts  on  the  properties  in  each  of  the 
expressions  below  refer  to  the  reactions  numbered  as  in  (2.40),  and  the  subscript  i  refers  to  the 
individual  species  in  the  mixture. 


R1  = 


X 


2m  i 

X  GiT'K'Pi 


R2  = 


_  S 


2mi 


Cie'e°'b  -  —C\ 
PD\  . 


-  —c\ 


PDi 


X  GiT^Pi 

R3  =  — 


m3 


pD\mm2 1 


C\C2  -  cie 


■^fr 


(2.41) 
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R4 = LaU + jd— _ 

m4  [  w*4m2/  Kc\mzm\ 

(ri)  +  W 

ms  l\msm\l  %  \ 


<PiT^p2e 


ig-^Difr 


\m$m\ 


m^mi 


These  expressions  for  the  non-dimensional  rates  of  each  reaction  have  been  written  without  the 
starred  notation  for  clarity.  Notice  that  the  subscripts  on  &d  and  po  refer  to  properties  associ¬ 
ated  with  each  reaction  and  its  Arrhenius  rate  expression.  Both  of  the  exchange  reactions 
contain  the  equilibrium  constant  Kc.  Finally,  m  for  each  reaction  retains  its  previous  definition 
as  the  average  mass  per  nucleus  of  the  reactants  normalized  by  the  atomic  mass  of  atomic 
hydrogen. 


The  source  term,  Wj ,  for  any  species  follows  upon  summing  the  contributions  from  all 
reaction  rate  expressions  affecting  that  species. 

Wx  =  *  mi(-2Rl  -  R3  -  R4  -  R5) 

W2  -  =  m2(-2R2  -  R3  +  R4  +  R5) 

R3-R4  +  R5>  (2.42) 

=  m4(Rl  +  R4) 

tf5  =  ^  =  m5(R2-R5) 


As  mentioned  previously,  Equation  (2.35)  eliminates  the  need  for  the  last  two  source 
expressions. 
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2.5  Rate  Constants 


Unless  specifically  noted,  all  calculations  completed  for  Part  III  use  rate  constants  and 
temperature  exponents  as  presented  in  Reference  30. 
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II  AN  ADAPTIVE  NUMERICAL 
METHOD 


3 .  Numerical  Integration  of  Governing  Equations 


Several  levels  of  embedding  removes  all  structure  from  a  uniform  initial  mesh.  Such 
unstructured  grids  do  not  fit  neatly  into  rectangular  matrices  and  demand  unstructured  data 
storage.  A  numerical  integrator  for  such  an  environment  should  operate  on  a  cell  by  cell  basis, 
dependent  only  upon  information  contained  within  a  single  cell.  This  effectively  rules  out  most 
implicit  schemes  and  focuses  attention  on  explicit  integration  schemes. 

This  need  for  a  compact  computational  stencil  resulted  in  the  selection  of  Ni’s  (29)  finite 
volume  formulation  of  the  classic  Lax-Wendroff,  explicit,  time  marching  scheme.  Pervaiz  (33) 
extended  this  scheme  to  include  chemical  source  terms  in  1987  and  Kallinderis  and  Baron  (16) 
recently  formulated  a  consistent  viscous  integration  extension.  The  algorithm  stores  the  state 
vector  at  cell  vertices  and  integrates  each  cell  independently,  requiring  only  cell-based 
information,  and  preserving  second-order  accuracy  at  computational  nodes. 

This  chapter  briefly  describes  inviscid  integration  on  2-D  cartesian  meshes  and  includes 
notes  concerning  chemical  source  terms,  smoothing,  and  boundary  conditions.  For  complete¬ 
ness  and  convenience,  Appendix  A  contains  the  additional  equations  for  2-D  non-orthogonal 
curvilinear  coordinates  (Appendix  A-l),  axisymmetric  non-orthogonal  curvilinear  coordinates 
(Appendix  A-2),  and  viscous  integration  in  non-orthogonal  2-D  curvilinear  coordinates 
(Appendix  A-3). 
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3.1  Integration  in  One  Dimension 


L 

R 

i  - . 

1  i 

i  + 1 

FIGURE  3.1 

A  one-dimensional  computational  domain. 


Consider  the  one-dimensionai  domain  of 
Figure  3.1.  The  1-D  governing  equation 
system  in  strong  conservation  form  is 

U,  =  -F,  +  W  ai) 


and  the  change  at  node  i  occurring  after  Ai  is 
SU/sU"+1  -U"  = 


U,  At  +  U„  ^  +  0(A/3) 


(3.2) 


Dropping  terms  of  third  order  and  higher  and  substituting  from  Equation  (3.1)  for  U,  provides 
the  second-order  total  change  at  node  i  based  on  contributions  at  the  n*  time  level. 


su,  =  A((W  -  F,r  +  ^(W0(W  -  F.HMW  -  FJU  I" 


(3.3) 


Equation  (3.3)  contains  the  Lax-Wendroff  step  and  evaluates  the  second  derivatives  of  the  state 
vector  using  the  flux  Jacobians  W u  and  Fy. 

U„  =  Wt/U<-(Ft/U<), 

The  first  and  second  terms  in  (3.3)  contain  the  first-  and  second-order  changes 
to  node  i  .  Ni's  primary  contribution  involved  reinterpreting  these  terms  based  on  cell- 
centered  values  of  the  state,  flux  and  source  vectors.  To  see  this,  recognize  that  the  source 
vector  at  /  is 

*  ^  4\ 


Here  W l  and  Wr  represent  cell  centered  values  of  W.  Using  central  differencing 


allows  reexpressing  the  first  parenthesis  of  (3.3)  in  terms  of  cell-centered  values. 

A/fW-F.l^AUt  +  aU,) 

Here  AUz,  and  AUl  are  the  cell  changes  defined  by: 

AUj^Wj^-IFi-F,,,)^. 


Similarly,  we  may  recast  the  second-order  term  of  (3.3)  using  these  cell  changes. 

W^W  -  Fjt)  - [¥u(W  -  Fx)]x }  =  f  {^(AU/,  +  AU/f)  -  [Ft,(AUL  +  AU R)l  J 

Defining  the  flux  change  and  source  change  vectors 

AF/,s  FyAUf, 

AWisWf/AUi 


allows  the  second-order  term  of  (3.3)  to  be  expressed  simply  as 


A/ 

2 


AW/,  +  AW*  AF/.-AF/} 
2  Ax 


(3.6) 


Combining  Equations  (3.5)  and  (3.6)  to  form  the  total  change  at  i  then  results  in: 

(5U,  =  i|AU  +  — AF+  ^Aw)  +  l{AU  -  ^AF+  ^Aw) 

2\  Ax  2  1L  2\  Ax  2  Jr  (3.7) 

This  form  clearly  shows  the  origin  of  contributions  to  node  i  from  both  adjacent  cells. 
This  property  makes  Ni's  scheme  very  attractive  for  use  with  unstructured  meshes.  Changes 
to  each  node  may  be  computed  separately  for  each  cell.  Then  all  contributions  to  any  node 
simply  may  be  summed  according  to  Equation  (3.7). 

More  precisely,  with  Equation  (3.7)  written  as 
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the  contributions  to  i  from  each  cell 

5uHAU+£AF+fH 

8U*,  =  (aU-^-AF+^Aw) 

'  1  Ax  2  )R  (3.9) 

may  be  computed  on  a  cell-by-cell  basis.  These  distribution  formulae  make  use  of  cell  center 
and  nodal  information  only,  fulfilling  the  stated  requirement  for  only  cell-based  information. 


3.2  Ni  Scheme  in  Two  Dimensions  with  Source  Terms 

In  two  spatial  dimensions,  the  integration  scheme  closely  resembles  the  one-dimensional 
formulation.  The  governing  equations  now  contain  a  second  flux  vector  as  shown  here. 

U«  =  -(FX  +  Gj,)  +  W  (3.10) 


An  analysis  comparable  to  that  of  the  last  section  leads  to  corresponding  nodal  change  formulae 
and  auxiliary  equations.  On  uniform  cartesian  meshes,  like  that  in  Figure  3.2,  the  change  at 
node  i  is 


8Uj 

ft 

FIGURE  3.2 

A  two-dimensional  computational  domain. 


AU^£AF^AG*  +  fAW* 

+  (aU»  -  |£aFj+  ^AG»  +  yAWg 
+  (AUc-^AFc-^AGc  +  fAWc 

+(auo  +  ^afo-0ago+|awo; 


(3.11) 


where  AFc  =  Ft/AUc 
AGc  —  G(/AUc 

AWc  =  Wt/AUc 


(3.12) 
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Rather  than  repeat  the  previous  analysis  to  derive  (3. 1 1),  we  consider  an  equivalent  finite 
volume  approach.  This  alternate  derivation  provides  greater  physical  insight  into  the  behavior 
of  the  scheme.  Here,  we  intend  to  examine  the  2-D  equations,  and  also  show  how  the  distribu¬ 
tion  formulae  assemble  the  governing  equations.  The  Lax-Wendroff  scheme  is  second  order 
accurate  in  time  and  space.  Ni's  formulation  (without  dissipation)  retains  this  order  of  accu¬ 
racy  at  the  nodes,  despite  the  fact  that  each  cell  is  computed  separately.  The  advantage  of  this 
property  cannot  be  overstressed  in  the  context  of  the  present  unstructured  solver. 


In  Equation  (3.9)  the  contribution  to  each  node  consists  of  primary  changes  from 
neighboring  cells  plus  or  minus  secondary  corrections  for  the  fluxes  through  cell  boundaries. 
In  terms  of  a  difference  stencil  at  any  node,  it  is  clear  that  each  cell  on  either  side  of  a  node 
performs  half  of  the  standard  central  difference  across  that  node  and,  when  added  together  at 
each  node,  these  changes  complete  the  stencil.  As  a  result,  in  the  steady  state,  AU  for  any 
particular  cell  may  be  non-zero,  provided  it  cancels  exactly  with  the  AU  contribution  from  its 
neighbor  at  the  nodes. 


In  keeping  with  a  cell-based  approach,  the  differential  governing  equation  (3.10)  should 
be  replaced  by  its  integral  form. 


(F*  +  G  y)dA 


I 


WdA 


The  divergence  theorem  changes  the  integral  of  flux  vectors  into  a  loop 
integral  over  the  cell's  perimeter. 


(Fdy  -  Gdx)  = 


WdAc 


Letting  U  and  W  take  on  cell-averaged  values  rids  the  far  left  and  right  terms  of  their  integral 
signs.  For  some  cell  C  on  the  mesh  shown  in  Figure  3.3, 
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(3.13) 


U,+ 


(F dy-  Gdx)  =  W 


Over  a  short  period  of  time  At  Equation  (3.13)  predicts  the  first  order  change  in  C. 


AUc  =  A/Wc- 


(F dy  -  Gdx) 


This  form  displays  the  nature  of  the  original  conservation 
equations.  The  first-order  change  in  any  cell  results  from  the 
average  contributions  of  sources  within  the  cell  and  the  flux 
through  its  boundaries.  Expanding  the  surface  integral  using 
the  midpoint  rule  gives 


FIGURE  3.3 

A  two-dimensional  computational  domain. 


AUc  =  AfWc  — 


At 

Ac 


+ 


m 

Vi  yi)  [—i-j. — 

(3.14) 


Close  examination  reveals  that  this  equation  reduces  exactly  to  the  "cell  changes"  defined 
after  Equation  (3.5).  However,  the  physics  exposed  by  the  conservation  law  of  (3.13)  offers 
strong  motivation  for  defining  AUc  in  such  a  way. 

With  this  in  mind,  return  to  the  differential  Equation  (3.10).  Substituting  the  Lax- 
Wendroff  step  from  (3.2)  into  (3.10)  results  in 

8Ui  =  aur  +  |(aw,-[AF,],  -[AG,],  r  (3  15) 

Here  AU);  and  AW);  are  cell  average  quantities  for  cells  A  •  D  (Fig.  3.3) 
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AUfc  =  +  AUfl  +  AUc  +  AU0) 

AW)/  =  ^(AW*  +  AWa  +  AWC  +  AWd) 


(3.16) 


and  the  derivatives  of  the  flux  Jacobians  are 

[AF/L  =  l/^gazMd.  +  AF_c-AFd\ 
21  Ax  A xl 

[AG/L  =  l/AGa-AG^  +  AGc  ~  AGp 

21  Ajc  Ax 


Back  substituting  Equations  (3.16)  and  (3.17)  into  (3.15)  results  in 


+ 


8U,  =  1 


+ 


(au^AF^AG^dw,) 

K-^AF»+^AG»  +  fAW»j 
K-^AFc-0AGc  +  fAWc) 
(au0  +  ^af0-^ag0+|aWd) 


which  is  exactly  Equation  (3.1 1). 


(3.17) 


(3.18) 


Again,  we  see  that  each  node  receives  a  contribution  from  surrounding  cells  to  complete 
the  differencing  and  preserve  second-order  accuracy  at  the  nodes. 

The  sketch  on  the  left  of  Figure  3.4  shows  a  schematic  of  the  changes  accumulated  at 
some  node  as  described  by  (3. 1 8). 


The  right  side  of  this  same  Figure  presents  the  same  information  from  a  cell-based  point 
of  view.  Each  cell  distributes  to  its  nodes  and  we  may  rewrite  Equation  (3.18)  from  this  cell- 
based  perspective,  giving  the  distribution  formulae  for  the  2-D  scheme. 
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FIGURE  3.4 

Schematic  of  changes  accumulated  at  the  Ith  node  in  a  two-dimensional  domain. 


Stic)/  =  $[aUc  -  £afc-  I^AGc  +  y&wc 

SUcI  -  J[aUc  +  ^AFc-^AGc  +  yAW C 

SUc)t  =  J{aUc  +  £-AFc+  0 AG  C  +  yAWc 

SUCU  if AUc  -  -~AFc+  ~AGC  +  ^AWC 

'  41  to  4y  2  j  (3.19) 

Equation  (3.19)  implicitly  assumes  identical  time  steps  for  each  of  the  cells.  However, 
for  cells  of  similar  size,  local  time  stepping  does  not  present  a  problem  in  steady  state 
calculations.  The  time-accurate  work  of  Pervaiz  (33)  discusses  these  points  in  some  detail. 

In  the  preceding  discussion  AF,  AG,  and  AW  involve  the  product  of  the  Jacobian 
matrices  and  changes  to  the  state  vector.  For  finite  rate  reacting  systems,  these  terms  differ 
from  the  original  perfect  gas  modeling  presented  by  Ni  (29).  Appendix  B  contains  these 
matrices  for  a  calorically  perfect  gas  (B-l),  and  the  multiply-reacting,  nonequilibrium  gas 
detailed  in  the  previous  section  (B-2).  For  gas  models  with  more  elaborate  modeling  of 
internal  energy  modes  and  energy  transfer  processes,  these  matrices  may  not  always  exist  in 
analytic  form. 
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3.3  Smoothing  Formulation 


Although  the  Euler  equations  do  not  include  any  dissipation,  the  Lax-Wendroff  discrete 
approximation  does.  The  scheme’s  truncation  error  introduces  fourth-order  terms  into  the 
modified  PDEs  capable  of  dissipating  weak  waves.  Figure  3.5  displays  this  property  through 
a  plot  of  the  amplification  factor  of  the  finite  difference  Equations  (2).  While  this  behavior  may 
provide  enough  damping  to  attenuate  weak  waves,  it  can  be  insufficient  for  the  steeper 
gradients  typically  found  in  hypersonic  calculations.  Additionally,  the  central  differenced 
algorithm  permits  odd-even  (sawtooth)  waves  to  remain  undetected  and  these  waves  must  be 
damped. 


Amplification  Factor  IGI 

FIGURE  3.5 

Amplification  factor  of  the  Lax-Wendroff  finite  difference  scheme.  (From  Ref.  2.) 


Smoothing  Formulation  in  One  or  Two  Spatial  Dimensions 

Referring  to  the  cells  in  Figure  3.6  and  taking  a  as  the  artificial  dissipation  coefficient,  the 
artificial  dissipation  at  node  i  is 

Di  -  f  (f|u)  *  f 

Or,  in  terms  of  the  average  value  of  the  property  u  in  each  cell 
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Dj  =  -&-(Ha  -  u„)  +-G-{Hb  -  un) 
A$  A«5 


FIGURE  3.6 

A  one-dimensional  computational  domain  smoothing 
formulation. 


This  formulation  extends  easily  into  two  dimensions.  Referring  to  Figure  3.7,  the  dissipation 
vector  at  ij  becomes 


(U.-UyMU.-Uy# 

+(Uc-u1'JMCo-iy/ 


(3.20) 


Defining 
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i 

D 
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C 
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B 

\ 

FIGURE  3.7 

A  two-dimensional  computational 
domain  smoothing  formulation. 


and  the  smoothing  contribution  to  the  i*  node  of  any  cell  C  as 

Dc=/i(Uc-U,) 

alters  the  two-dimensional  net  change  to  i  (Eq.  3.1 1)  as  follows. 

'  ( AU*  +  £&Fa  +  ^AGa  +  y  A  W A  +  Da 

+  (AUn  -  |^AFfl+  ~AGB  +  yAWfl  +  DB 
+  (aUc  -  |^-AF  c-  ^-AGc  +  y  AW  c  +  Dc 
+  (AUo  +  ^AF d-  |^AGp  +  yAWp  +  Dp 


8Uf  =  l 


(3.21) 


43 


The  shock  fit  approach  removes  the  bow  shock  from  the  domain  "interior,"  and  virtually  the 
entire  flowfield  is  smooth  by  comparison.  Nevertheless,  the  "smooth"  region  is  still  strongly 
nonlinear  and  requires  stabilization  throughout.  Placing  the  discontinuity  at  the  boundary  of 
the  computational  domain  effectively  removes  the  need  for  alternate  smoothing  levels  near  and 
away  from  discontinuities. 

Behavior  of  Dissipation  in  Goveming^quations 

Criticism  of  smoothing  operators  stems  from  the  fact  that  they  modify  the  original  set  of 
equations.  To  assess  how  the  second  difference  smoothing  operator  alters  the  actual  system, 
examine  the  change  at  any  node  i  in  a  one-dimensional  domain. 

AUj  +  I^AFj  +  yAW„  +  gy(--y-"‘:) 

4Uc_£AFc  +  |aWc+o|(i!i±J_^) 

FIGURE  3.8 

A  one -dimensional  computational  domain  (3.22) 

for  examining  smoothing  contributions. 


13 

c 

_ 4 

1 - ►, 

i  •  I 


SU,  =  1 


This  is  Equation  (3.21)  for  a  single  dimension.  The  first  bracket  contains  the  change  to  i  from 
B  and  the  second  contributes  the  change  from  C.  Decomposing  this  equation  further  results  in 
a  modified  form  of  the  Lax-Wendroff  step  (Equation  3.3). 

8U,-  =  Ar(W  -  FJ  +  ^{Wu(W  -  F,)-[Fu(W  -  ¥x)l }  +crAf  (U)** 


Equations  (3.3)  an  (3.23)  differ  only  by  the  smoothing  term  on  the  extreme  right  of  the  latter 
expression,  forcing  us  to  examine  the  magnitude  of  the  error  introduced  by  this  term. 

(<t82u)a/ 


8U;  = 


A/2 

U  tAt  +  U  tt~\ 


(3.24) 


Here  operator  notation  has  been  used  to  highlight  the  second  difference  smoothing.  Since  the 
smoothing  adds  a  first-order  error  to  the  change  at  i,  cr  must  be  kept  small  enough  to  avoid 
changing  the  value  of  the  state  vector  appreciably.  Fortunately,  in  smooth  regions,  the  second 
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difference  of  U  tends  to  vanish  and  the  scheme  retains  its  accuracy  even  for  moderate  values  of 
<r.  Equation  (3.24)  also  predicts  the  algorithm  will  contain  the  largest  errors  where  flow 
gradients  change  most  rapidly.  In  practice  a  varies  from  0.1  to  1.0  for  hypersonic  calculations 
and  from  about  0.02  to  0. 1  for  transonic  and  subsonic  problems. 


The  form  of  Equation  (3.22)  demonstrates  that  the  overall  smoothing  contribution  to  any 
node  results  from  taking  the  difference  of  two  first  differences.  That  is 

(ui+i  -  Ui)-(ui  -  U(_i)  =  Aui+y2 -  Au/_y2  =  m~i  -  2m  +  ui+i 
Each  cell  contributes  one  first  difference,  forming  the  second  difference  when  the  changes  from 


both  cells  combine  at  any  node. 


FIGURE  3.9 

Variation  of  property  u  showing  degeneracy  of 
smoothing  operator  near  boundaries. 


Consider  the  special  case  of  a  linearly  varying 
property  u  across  nodes  ij,  to  k  at  the  boundary. 
Since  the  second  difference  measures  curvature,  it 
should  return  a  zero  at  any  of  these  nodes.  In 
particular,  the  smoothing  contributions  from  cells  A 
and  B  to  nodes  j  and  k  are: 


at  the  interior  node./;  Dj  =  DA)j  +  DB)j  =  f*  +  =  0 

(3.25) 

at  the  boundary  node  k\  0*  =  0a)*  +  O  =  (^  +  Oj^*O 

Since  the  operator  detects  curvature  and  the  line  is  straight,  zero  contributions  are  expected  at 
both  nodes.  However,  only  cell  B  contributes  to  the  smoothing  stencil  at  k,  and  a  first 
difference  operator  results.  While  the  operator  behaves  correctly  on  interior  points,  it 
degenerates  along  boundaries.  Notice  that  the  first  difference  operator  appears  as  a  direct  result 
of  the  boundary  and  is  not  related  to  the  property  variation. 
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Generalizing,  the  modified  governing  equations  take  the  following  forms: 


On  interior  nodes 
Along  £  =  const,  boundary 
Along  tj  =  const,  boundary 


U«  +  F$  +  Gtj  =  W  +  +  8^u) 

U„  +  F|  +  Gn  =  W  +  +  A^u) 

U„  +  F4  +  G,,  *  W  +  o(a4U  +  S^u) 


(3.26) 


Obviously,  the  first  differences  along  the  boundaries  will  create  severe  errors  in  flows  with 
large  spatial  gradients.  However,  since  the  second  difference  operator  measures  curvature,  it 
will  remain  small  in  such  fields  -  provided  the  gradients  themselves  do  not  change  rapidly. 
Additionally,  while  the  second  difference  operator  is  purely  dissipative,  the  first  difference 
behaves  convectively,  significantly  changing  the  nature  of  the  smoothing  terms. 

Not  smoothing  in  a  direction  normal  to  the  boundaries  avoids  this  error.  Dropping  the 
first  difference  terms  on  the  right  side  of  Equation  (3.26)  results  in: 

On  interior  nodes  U„  +  F$  +  G,,  *  W  +  o(8$U  +  8ju) 

Along  £  =  Const,  boundary  U„  +  F^  +  Gn  =  W  +  o^U )  (3.27) 

Along  rj  =  Const,  boundary  U„  +  F$  +  G  n  =  W  +  ) 

Stagnation  Enthalpy  Smoothing 

Writing  the  equations  as  in  (3.27)  shows  that  the  smoothing  terms  behave  like  source 
terms  in  the  governing  equations.  To  minimize  changing  the  physics  of  the  problem,  we 
choose  <xto  keep  such  "sources"  small. 


Without  smoothing,  the  steady-state,  inviscid  energy  equation  holds  stagnation  enthalpy  con¬ 
stant.  In  one  dimension 


Be 

Bi  +  pU 


=  o 


(3.28) 
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In  the  steady  state  this  becomes 

pu<*M  =  0 

(3.29) 

Anv  additional  terms  violates  this  adiabatic  statement.  For  a  nonequilibrium  system, 
enthalpy  is  representative  of  many  modes  and  an  altered  stagnation  enthalpy  may  lead  to 
unrealistic  solutions.  For  example,  if  a  nonequilibrium  mode  absorbs  only  1%  of  the  total 
enthalpy,  and  the  error  in  stagnation  enthalpy  is  1%,  that  mode's  behavior  may  be  totally 
incorrect.  Although  the  net  effect  will  be  small  (only  0.01  h0),  the  outcome  may  mask  the  true 
physics  of  the  problem.  Since  hypersonic  flows  are  characteristically  high  enthalpy  flows, 
these  observations  can  have  serious  implications. 

Errors  in  stagnation  enthalpy  have  other,  more  subtle,  effects.  Consider  the  schematic  in 
Figure  3.10.  For  a  stagnated,  frozen,  perfect  gas,  all  enthalpy  must  reside  in  random  particle 
motion,  and  constant  stagnation  enthalpy  implies  constant  stagnation  temperature.  When  the 
flow  enters  the  shock  layer  at  point  b,  it  very  nearly  stagnates  behind  the  normal  shock.  The 
post-shock  temperature  at  b'  jumps  to  within  a  few  percent  of  T0.  In  slowing  isentropically  to 
the  stagnation  point  c,  T  increases  very  little.  However,  this  small  variation  is  of  primary 
interest  within  the  shock  layer  and  even  a  small  error  in  total  enthalpy  can  appear  large  by 
comparison. 


FIGURE  3.10 

Stagnation  streamline  temperature  behavior  in  a  frozen  flow. 
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The  one-dimensional  energy  equation  with  artificial  smoothing  terms 


*+*£!*!!>.  0(A) 

3<  as 


(3.30) 


no  longer  implies  constant  hQ.  Instead,  the  balance  forces  total  enthalpy  to  vary  in  space  in 
accord  with  the  non-linearity  of  the  total  internal  energy.  The  term  is  small  but  the  preceding 
discussion  indicates  the  possible  consequences.  In  practice,  a  a  level  of  0.1  has  a  noticeable 
effect  on  h0  in  the  domain  downstream  of  the  shock. 


To  prevent  this,  a  second  smoothing  term  adds  a  correction  to  the  internal  energy  at 
i  to  counteract  the  error  produced  in  (3.30).  This  correction  forces  h0i  toward  h0 tending  to 
hold  ho  constant  throughout  the  domain.  Specifically, 


P  y  _  ^K,{ho.  ~  h0i)ej 
1  *  ho. 


(3.31) 


Here  the  subscript  i  denotes  local  quantities,  Pho  is  a  stagnation  enthalpy  smoothing  coefficient 
and  the  ratio  eilho-  rescales  the  term  for  use  in  the  energy  equation.  If  h0  =  const  -  h0oo  this 
correction  vanishes.  In  other  words,  the  term  contributes  only  if  the  standard  second  differ¬ 
ence  smoothing  violates  the  physics  expressed  by  the  energy  equation.  Adding  this  term  to  the 
energy  equation  (only)  results  in  the  following  change  in  total  internal  energy  at  any  node  i. 


'  (aU„  +  £aF„  +  I^AG*  +  +  D„  +  04.) 

+  (aUs  -  f^AF„+  ^AG„  +  yAWs+  Dfl  +  04) 

+  (aUc  -  ^AFC-  |^AGc  +  yAWc  +  DC  +  DjJiJ 
+  (aUo  +  ^AFo-  ^AG0  +  ^AWd+Dd  +  04) 


In  practice  a  value  of  pho  two  to  three  times  larger  than  the  overall  smoothing  coefficient  p 
suppresses  total  enthalpy  variations  to  0.1%  or  less.  Use  of  such  a  large  coefficient  seems 
reasonable  since  the  term  tends  to  restore  the  physics  of  the  original  governing  equations. 


As  a  final  comment  on  stagnation  enthalpy  smoothing,  note  that  smoothing  terms  affect 
all  of  the  conservation  equations.  At  any  point  in  the  field  they  effectively  create  mass  sources, 
momentum  sources,  etc..  Reference  7  shows  that  the  second  difference  operator  is  globally 
conservative  for  constant  smoothing  coefficients.  Since  the  other  conservation  statements 
(mass,  momentum,  species)  are  integral  properties,  the  local  sources  and  sinks  globally  cancel 
due  to  the  conservative  nature  of  this  smoothing  formulation.  On  the  other  hand,  stagnation 
enthalpy  is  a  point  property,  and  while  overall  eneigy  may  be  conserved,  any  source  term  will 
be  incorrect  at  a  particular  point 


3.4  Numerical  and  Physical  Boundary  Conditions 

A  typical  hypersonic  shock  layer  maps  from  physical  to  computational  space  as  shown 
below  (Figure  3. 1 1).  All  flow  in  the  domain  passes  through  the  shock  (ad)  and,  in  the  steady 
state,  must  be  balanced  by  the  out  flow  across  the  shock  layer  ( cd ).  When  the  gas  is  assumed 
inviscid,  it  slips  tangentially  over  the  body's  surface. 


Plane 

FIGURE  3.11 


Schematic  of  physical  and  computational  space  for  blunt  body  computations. 
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For  symmetric  two-dimensional  problems,  the  natural  symmetry  plane  allows  calculation 
of  only  half  the  domain.  In  the  axisymmetric  case,  the  stagnation  streamline  becomes  the  axis 
of  symmetry,  again  permitting  this  simplification. 


Edge  Cells  in  Finite  Volume  Schemes 

Figure  3.1 1  contains  four  types  of  boundaries  surrounding  the  computational  mesh,  and 
all  boundary  nodes  have  physical  cells  on  only  one  side.  For  cell  vertex  based  schemes  this 
results  in  incomplete  differencing  across  such  nodes. 

Consider  the  boundary  cells  A  and  B  in  the 
sketch  shown  to  the  left  (Figure  3.12).  The 
nodes  i,  j,  and  k  receive  corrections  8 A  and 
8B  from  one  side  only.  In  other  words, 
boundary  nodes,  such  as  these,  receive  only 
two  corrections,  while  interior  nodes  receive 
four. 


\  / 

A 

/  \ 

\  7 

A' 

\  y  >  >  n  \  \  y 

x  7 

B' 

FIGURE  3.12 

Edge  cells  in  a  node-based  finite-volume  scheme. 


From  a  finite  difference  standpoint,  the  spatial  derivatives  make  use  of  a  central  difference 
operator  at  all  computational  nodes.  Normal  to  boundaries,  this  operator  becomes  one-sided 
and  requires  special  treatment. 

To  prevent  this  degeneracy,  imagine  a  second  set  of  cells  A'  and  B‘  outside  the  physical 
boundary.  These  fictitious  cells  complete  the  central  difference  taken  normal  to  the  boundary 
and  permit  upgrading  of  these  nodes.  In  this  figure,  cells  A'  and  B'  contribute  a  change  S/T 
and  SB'  to  j  providing  it  with  a  total  of  four  changes  -  just  like  an  interior  node.  In  the  event 
no  better  estimate  exists,  setting  &4'  and  6fl'  equal  to  &4  and  SB  results  in  first  order  accuracy 
(since  this  is  equivalent  to  a  forward  or  backward  difference). 


Implementing  this  general  treatment  requires  simply  doubling  the  changes  to  each 
boundary  node.  Each  of  the  boundaries  shown  in  Figure  3.1 1  presents  a  physical  situation 
which  permits  modification  of  this  procedure  to  improve  the  accuracy  of  this  procedure  based 
on  knowledge  of  the  physical  boundary. 

Inflow  /  Shock  Boundary 

The  bow  shock  forms  the  inflow  boundary  of  the  domain.  This  confines  the  domain  to 
points  within  the  shock  layer  and  avoids  introducing  undisturbed  cells.  Chapter  S  describes  the 
shock  fitting  procedure  in  some  detail. 

Symmetry  Plane  Boundary 

The  symmetry  plane  presents  precisely  the  situation  shown  in  3.13.  The  imaginary  cells 
A'  and  B'  are  simply  mirror  images  of  A  and  B.  The  boundary  conditions  for  such  nodes  is 
clear.  Using  the  notation  of  Figure  3.12,  the  change  at  node;  is: 

8U;  =  28U*)y-  +  28Ufl)y  (3.34) 

Since  it  is  also  a  streamline  of  the  flow,  the  normal  velocity  vanishes  identically.  This 
recognizes  no  flow  across  the  boundary  and  prevents  the  growth  of  spurious  error  and 
sawtooth  oscillations. 
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Outflow  Boundary 


Schematic  of  outflow  boundary  in  bhmt  body  flow  showing  location  of  the  sonic  line. 


Doubling  the  changes  for  the  exit  plane  results  in  only  first  order  accuracy  across  the 
downstream  boundary.  However,  for  most  configurations,  the  supersonic  flow  across  this 
plane  does  not  permit  the  error  to  propagate  upstream  and  corrupt  the  rest  of  the  domain. 


Solid  Wall  Boundary 

Inviscid  flow  slips  tangentially  along  the  wall  following  the  local  surface  inclination  and 
the  body  surface  streamline.  For  low  enthalpy  flows,  this  boundary  condition  is  easily  imple¬ 
mented.  At  each  time  step,  one  computes  the  total  velocity  at  a  body  node,  discarding  the  nor¬ 
mal  components.  Alternatively,  when  integrating  the  cells  along  the  solid  wall  boundary,  one 
may  enforce  a  no-flux  condition  along  the  body  surface.  This  prevents  the  face  from  con¬ 
tributing  to  the  flux  integral.  These  simple  treatments  usually  perform  adequately.  However, 
the  high-enthalpy  flows  presendy  under  consideration  require  greater  care. 


In  these  flows  •  especially  near  the  stagnation  region  •  these  treatments  may  lead  to  strong 
transients  during  convergence,  or  extreme  sensitivity  to  initial  conditions.  In  reacting  flows,  or 
gas  models  with  multiple  internal  energy  modes,  any  error  in  stagnation  enthalpy  may  lead  to 
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non-physical  solutions.  Clearly,  discarding  the  normal  velocity  component  throws  away 
energy.  To  avoid  such  enthalpy  errors,  it  is  better  to  reorient  the  velocity  vector  at  each  node. 


FIGURE  3.14 

Flow  tangency  condition  for  inviscid  simulations  showing  reorientation  of  the  total  velocity 
vector. 


The  velocity  vector  V  consists  of  the  current  (time  n)  values  of  velocity  plus  changes  to  this 
vector  8V.  After  first  computing  this  quantity  at  some  surface  node  i. 


\Vi\*W+W 

with  x  and  y  components 

u}2=  u ?  +  28u, 
Vi2=  vf  +  28v,- 


(3.35) 


rotate  this  vector  to  the  local  surface  slope  inclination  angle  to  preserve  translational  energy. 


«i 

Vi* 


7*1*1 +  7*2*1 

KK  ‘d|dS 

7*2*1+ 7*2*1 


(3.36) 


Finally,  rescale  the  resulting  vector  to  the  same  overall  magnitude  of  the  original  velocity  vec¬ 
tor. 


(3.37) 
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At  first  glance,  forcing  the  velocity  vectors  tangential  at  each  body  node  may  appear  equivalent 
to  enforcing  the  no-flux  condition  on  the  surface.  The  nodal  velocity  will  convect  members  of 
the  state  vector  parallel  to  the  surface  preventing  these  fluxes  from  contributing  to  the  integral 
performed  upon  each  cell. 


Certainly  this  is  the  case  for  both  linear  and  parabolic  bodies,  but,  upon  more  generally 
shaped  bodies,  flow  tangency  at  surface  nodes  may  not  entirely  prevent  contributions  to  the 
flux  integral.  The  cell's  wall  may  not  necessarily  be  at  the  same  (algebraic  mean)  slope  of  the 
comer  nodes.  While  this  is  a  higher  order  error,  experience  suggests  that  not  enforcing  a  no¬ 
flux  boundary  may  degrade  the  initial  convergence  properties  of  the  scheme. 


Examine  the  integral  expression  for  the  changes  to  cell  C. 


MJ)c  =  ArWc- 


F dy-  Gdx 


(3.38) 


The  contribution  to  the  loop  integral  over  the  south  face  is: 


FIGURE  3.15 

No  contribution  to  flux  integral  from  surface  face. 
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Re-arranging  yields 


FAy,  -  GAx,  = 


p 

0 

pu 

{uAy,  -  vAx,)  + 

pAy, 

pv 

-p  Ajc, 

ph. 

0 

.  A  . 

0  . 

(3.40) 


But,  no-flux  implies  n  •  V  =  0  through  the  southern  face. 

n.y.-fe4fr-vA**)sQ 


So, 


uAys  -  vAx,  =  0 


(3.41) 


Equation  (3.41)  forces  the  entire  first  term  on  the  right  of  (3.40)  to  zero  leaving  only  the 
pressure  forces  in  the  second  and  third  elements  of  the  flux  integral  to  contribute  to  the  integral. 


FAy,  -  GAx,  = 


o 

pAy, 

-P&x, 

0 

0 


(3.42) 


Applying  Equation  (3.42)  during  the  integration  step,  reorienting  the  velocity  vector  by  (3.37), 
and  doubling  the  changes  to  the  other  members  of  the  state  vector  by  (3.37)  insures  correct 
boundary  treatments  and  avoids  most  problems  with  initial  transients. 
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4.  Adaptation  and  Unstructured  Meshes 


The  discrete  equations  only  approximate  the  governing  equations  of  inviscid,  reacting, 
supersonic  flow.  The  Lax-Wendroff  integration  scheme  neglects  terms  of  third  order  and 
above  in  both  space  and  time.  As  a  result,  it  may  become  necessary  to  increase  the  quality  of 
any  given  solution  by  increasing  the  grid's  spatial  resolution.  However,  since  smooth  and 
slowly  varying  regions  of  the  flow  field  contain  little  third  order  activity,  refinement  in  those 
areas  adds  little  to  the  global  accuracy  of  the  solution.  In  more  steeply  varying  regions, 
however,  the  linear  and  second  order  terms  in  the  FDE’s  do  not  permit  adequate  flexibility  to 
resolve  the  physics  of  compressible  flows.  As  discussed  previously,  this  situation  results  in 
the  largest  numerical  errors  where  the  flow  gradients  vary  the  fastest.  Unfortunately,  it  is 
precisely  such  physical  structures  which  often  dominate  a  flow’s  behavior. 

Grid  adaptation  attempts  to  resolve  this  disparity  by  either  re-distributing  or  embedding 
nodes  to  improve  resolution  of  flow  features.  References  9,  35,  and  7  contain  discussions 
comparing  the  relative  merits  of  both  approaches. 

Originally  developed  by  Dannenhoffer  and  Baron  (6),  the  technique  described  here 
refines  a  solution  by  embedding  finer  grids  in  response  to  flow  features  detected  on  coai  ■  er 
grids.  While  this  method  reliably  detects  these  features,  its  real  strength  is  the  ability  to  embed 
several  grid  levels  in  arbitrarily  shaped  regions  during  evolution  of  a  solution. 
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In  a  hypersonic  shock  layer,  one  expects  large  species  gradients  after  a  chemically  excited 
flow  passes  through  a  strong  bow  shock.  In  order  to  resolve  these  gradients,  the  chemical 
relaxation  length  cannot  be  much  smaller  than  about  one  percent  of  a  local  cell  dimension.  This 
effect  also  was  documented  by  Park  (30)  who  described  it  in  terms  of  the  mesh  Damkohler 
number.  Such  near  equilibrium  flows  present  significant  computational  difficulties  due  to  the 
small  length  scales.  Although  it  is  tempting  to  consider  identical  equilibrium  modeling,  the 
rapidly  varying  chemical  relaxation  length  demands  a  nonequilibrium  calculation  as  the  flow 
continues  through  the  shock  layer.  Dissociation,  for  example,  absorbs  much  of  the  flow's  en¬ 
ergy  and  may  change  the  chemical  length  scale  by  several  orders  of  magnitude.  In  addition  to 
providing  greater  resolution  of  the  relaxation  zone,  grid  adaptation  helps  to  resolve  behavior 
near  internal  shocks,  expansion  fans,  sliplines  etc.. 

Note  that  the  bow  shock  forms  the  inflow  boundary  and  the  domain  contains  no 
freestream  cells.  All  computational  cells  lie  in  "disturbed"  regions  and  contribute  to  the 
solution.  Thus,  shock  fitting  performs  a  somewhat  similar  function  as  adaptation.  Moreover, 
since  the  blunt  body  domain  is  narrowest  near  the  stagnation  point,  the  shock  layer  shape 
automatically  clusters  nodes  near  the  nose  during  the  initial  grid  generation  phase.  As  a  result 
we  expect  a  reasonably  applicable  coarse  mesh,  and  the  adaptation  serves  to  tailor  the  grid  by 
further  refining  structures  within  the  shock  layer. 


4.1  General  Procedure 

Adaptation  greatly  increases  the  quality  of  a  given  solution  with  modest  additional 
computational  effort.  The  results  published  by  Baron,  Dannenhoffer,  Kallinderis,  Pervaiz  and 
Shapiro  (6,  17,  33,  and  34)  demonstrate  this  quite  clearly.  Often,  equally  resolved  solutions 
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without  embedding  must  use  globally  fine  grids  and  are  too  expensive  to  compute,  even  to 
benchmark  a  particular  solution. 

This  added  efficiency  comes  at  the  expense  of  coding  effort.  To  some  degree,  the 
perceived  complexity  of  coding  an  unstructured,  adaptive  domain  prevents  the  widespread  use 
of  adaptive  techniques.  Here  we  describe  a  simple  (but  effective)  adaptive  algorithm. 
Additionally,  the  final  Section  (4.5)  stresses  the  ease  of  implementing  adaptation  routines  in  a 
reasonably  efficient  manner. 

The  basic  algorithm  contains  three  steps: 

/.  Detection:  Examination  of  the  fiowfield  in  search  of  flow  features. 

ii.  Division:  Dividing  large  cells  into  smaller  cells  to  create  new  cells  and 

nodes. 

iii.  Pointer  Updating:  Absorption  of  new  cells  and  nodes  into  the  existing  data 

structure 

While  the  procedure  described  here  is  similar  to  that  presented  by  Dannenhoffer  (6)  and  the 
spatial  adaptation  by  Pervaiz  (33),  the  implementation  is  believed  to  be  less  complex  than 
either.  As  a  result  it  lacks  some  of  the  subtler  features  of  these  previous  works.  Nevertheless, 
it  works  well  in  practice. 

4.2  Detection  of  Flow  Features 

Detecting  features  in  a  fiowfield  typically  requires  consideration  of  threshold  values, 
differences  in  computational  space,  and  independent  detection  parameters.  The  suggested 
scheme  examines  one  or  more  independent  properties  throughout  the  domain  and  evaluates 
either  first  or  second  differences.  Then,  after  normalizing  by  a  convenient  rule,  the  algorithm 
tags  cells  containing  differences  above  a  certain  threshold  for  refinement. 
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Feature  detection  is  not  a  unique  procedure,  and  while  it  is  relatively  simple  to  design  a 
system  capable  of  trapping  some  given  feature,  designing  a  detector  to  trap  ull  interesting  or 
important  features  is  a  much  more  difficult  task.  It  is  difficult  to  precisely  define  "interesting  or 
important."  More  importantly,  most  flow  features  do  not  have  clearly  defined  "edges,"  making 
threshold  selection  somewhat  arbitrary. 

Figure  4.1  indicates  some  of  the  expected  flowfield  structures.  These  features  result 
from  both  fluid  dynamic  and  thermo-chemical  phenomena  and  all  exist  to  varying  degrees  in 
most  hypersonic  flows.  The  inviscid  Euler  equations  predict  the  entropy  layer  and  expansion 
fan  in  supersonic  blunt  body  flows  for  gases  with  and  without  internal  degrees  of  freedom. 
The  nonequilibrium  features  arise  as  a  direct  result  of  a  real  gas  incorporating  into  the 
governing  equations. 


FIGURE  4.1 

Some  flow  feature*  in  hypertonic,  real -gas,  blunt  body  flows. 
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In  order  to  separate  the  gas  dynamic  and  chemical  features,  we  search  for  flow  variables 
which  respond  independently  to  fluid  and  nonequilibrium  chemical  structures.  As  an  element 
of  gas  flows  through  an  expansion  fan  e.g.,  p,  T  and  p  all  decrease  monotonically.  The 
smooth  behavior  implies  small  second  differences,  indicating  the  benefits  of  first  differences 
for  detection.  First  differences  furnish  a  slope  per  cell.  Over  a  typical  blunt  wedge,  p,  T  and  p 
tend  to  be  constant  in  both  the  stagnation  region  and  along  the  flat  body  surface  downstream  of 
the  expansion  fan,  and  differences  for  any  one  of  the  variables  will  correctly  identify  the 
feature.  T  and  p  respond  much  the  same  to  an  entropy  layer  and  either  of  those  properties 
would  identify  this  feature.  Reference  7  lists  several  gas  dynamic  features  and  evaluates 
criteria  for  locating  specific  structures. 

For  the  features  in  Figure  4. 1  the  choice  of  a  best  parameter  to  identify  the  fluid  dynamic 
structures  is  not  crucial.  Density  was  chosen  since  it  is  a  state  variable  and  mass  must  be 
conserved.  Moreover,  since  it  is  a  primitive  variable,  integration  yields  a  direct  value  without 
auxiliary  equations.  As  a  result  it  is  less  prone  to  spurious  errors  than  are  calculated  quantities 
like  pressure  or  temperature. 

The  choice  of  nonequilibrium  parameters  depends  largely  upon  the  internal  modes 
included  in  the  gas  model.  In  general,  these  modes  act  independently  and  on  different  physical 
scales,  and  each  requires  its  own  independent  parameter.  For  dissociating  and  vibrating  gas 
models,  both  species  concentrations  and  vibrational  temperature  should  be  considered  as 
independent  parameters.  For  cases  where  all  nonequilibrium  modes  are  chemical,  the  reaction 
rates  associated  with  production  of  different  species  vary  greatly  and  result  in  different 
reactions  occurring  in  various  parts  of  the  domain.  For  such  flows,  differences  in  all 
independent  species  indicate  chemical  behavior  along  the  various  reaction  paths. 


Of  course,  the  most  important  criteria  that  a  flowfield  indicator  must  satisfy  is 
independence.  Since  pressure,  temperature,  and  density  all  respond  to  an  expansion  fan,  any 
one  of  them  would  identify  the  feature.  Interrogating  more  than  one  is  unnecessarily 
expensive.  However,  through  that  same  expansion  fan,  the  chemistry  will  often  tend  to  freeze 
out  and  display  little  variation.  Since  it  responds  to  an  entirely  different  class  of  flow  feature, 
species  concentration  behaves  as  an  independent  parameter. 

Since  cells  spanned  by  large  differences  (in  computational  space)  divide  as  the  adaptive 
process  continues,  the  procedure  drives  the  solution  toward  uniform  differences  in 
computational  space.  Embedding  smooths  out  variations  in  computational  space.  The  first  and 
second  order  terms  preserved  by  the  Lax-Wendroff  differencing  permit  tracking  of  the 
governing  equations  up  to  the  second  derivatives.  This  smoother  domain,  with  its  smaller 
higher  order  derivatives  lends  itself  to  more  accurate  computation  of  the  state  vector. 


The  Mechanics  of  Detection 


FIGURE  4.2 
General  undivided  cell. 


For  a  general  cell  considered  by  a  feature  detector,  comer 
nodes  2, 4,  6  and  8  always  exist.  The  midface  nodes  3,  5, 
7,  and  9  may  exist  if  any  adjacent  cells  are  already  divided. 


The  difference  in  the  £  direction  is 


where  p  is  density  and  Cj  is  a  nonequilibrium  parameter  (i.e.  mass  fraction  of  the  1th  species). 
Similarly  the  difference  in  77  is 


8»jP 


A  rational  definition  gives  the  overall  change  in  any  cell;  e.g. 


(4.1) 


Ape  V (8{pf  +  (5 vpf 
Ac,  s  V(8$c/f  +  (Snc,j2 

These  definitions  could  easily  be  biased  toward  £  or  rj  by  altering  the  exponents  under  the 
radicals. 

The  magnitudes  of  the  differences  are  important  primarily  for  comparison  with  other  dif¬ 
ferences,  and  should  be  normalized.  Popular  choices  for  reference  levels  include  the  statistical 
mean,  mode,  or  maximum  value.  Choosing  the  maximum  value,  A p  max  and  Ac/max, 

conveniently  bounds  the  variation  of  flowfield  indices  between  zero  and  one.  Obviously,  the 
reference  choice  is  not  unique  and  most  parameters  work  equally  well. 

Adaptation  Maps  and  Thresholds 

After  computing  the  refinement 
parameters  for  each  cell,  it  is  helpful  to 
prepare  a  map  as  shown  (Fig.  4.3).  (The 
magnitudes  implied  by  the  positive  branch  of 
the  radical  in  Equation  (4.1)  actually  collapse 
the  map  to  the  1st  quadrant.)  This  example 
shows  a  map  for  two  parameters.  In  general, 
such  maps  contain  as  many  dimensions  as 
there  are  adaptation  parameters. 

Cells  containing  little  or  no  variation  will  cluster  near  the  origin,  while  cells  within  or  near  fea¬ 
tures  plot  lurther  from  the  origin.  Since  the  two  parameters  should  be  independent,  we  antici¬ 
pate  little  activity  along  the  45°  diagonal  Ac/  =  A p.  Clustering  along  such  a  line  would  indicate 
a  feature  being  tracked  by  both  parameters,  and  thus  parameters  which  do  not  respond 
independently. 

The  adaptation  parameter  measures  the  distance  from  the  origin  to  any  point  on  the  map, 


FIGURE  4 J 

An  adaptation  map  based  cm  two  parameters. 
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(4.2) 


and  the  scheme  marks  a  cell  for  division  whenever  AC  is  greater  than  some  threshold,  Ta^. 


For  n  adaptation  parameters,  AC  becomes  a  radius  in  n  dimensional  space,  and  Eq.(4.2) 
generalizes  to 


AC  =J{Apf  +  Y[Ac,f 


Again,  adjusting  the  exponents  under  the  radical  weights  any  parameter  accordingly. 

Figure  4.4  presents  an  example  of  an  actual  adaptation  map  computed  for  a  relatively  low 
enthalpy  Mach  5  flow  with  a  simple  dissociation  reaction  (A2  +  M  <=>  2A  +  M)  and  illustrates 
many  of  the  features  discussed  above. 


FIGURE  4.4 

Adaptation  map  for  Mach  S  flow  ihowing  polarization. 


Since  the  plot  exhibits  considerable  scatter,  it  is  fair  to  say  that  differences  in  various 
r  . -ions  remained  relatively  benign  with  respect  to  one  another.  More  precisely,  the  chemical 
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and  density  gradients  are  distributed  relatively  evenly  throughout  the  field  and  not  confined  to  a 
few  isolated  regions.  Since  nearly  all  cells  have  some  species  and  density  variation,  it  is  clear 
that  the  reactions  in  the  domain  did  not  occur  over  scales  very  disparate  from  the  physical  scale 
of  the  body.  For  example,  the  dissociative  activity  did  not  occur  solely  adjacent  to  the  shock. 

The  polarization  of  these  data  demonstrate  that  atom  mass  fraction  and  density  are 
responding  independently  as  intended.  This  polarization  cannot  be  overstressed  and 
demonstrates  that  each  of  the  adaptive  indices  responds  (independently)  to  its  own  type  of  flow 
features. 


FIGURE  4.5 

Behavior  of  a  general  adaptation  map. 


Figure  4.5  examines  this  behavior  in  more 
detail.  Here  a  mock  adaptation  map  shows 
ail  cells  having  AC  2  0.75  tagged  for  divi¬ 
sion.  Cells  in  region  A  demonstrate  such  ex¬ 
treme  nonequilibrium  that  they  will  be 
adapted  regardless  of  their  density  behavior. 
Typically,  cells  in  this  area  lie  in  rapidly  re¬ 
laxing  flow  regions.  Conversely,  cells  in  D 
show  large  density  differences,  but  little 


chemical  activity.  As  the  flow  expands  out  of  a  stagnation  region  and  around  a  body,  the 
density  drops  rapidly,  making  it  difficult  for  various  reactants  to  locate  collision  partners  and 
effectively  turning  off  the  chemistry.  Indicative  of  rapid  variation  in  both  parameters,  region  B 
contains  few  cells  in  practice.  However,  if  the  fluid  dynamical  parameters  depend  strongly  on 
the  degree  of  nonequilibrium,  cells  may  appear  here.  This  often  occurs  in  gases  such  as  ni¬ 
trogen  where  the  dissociation  energy  may  comprise  a  significant  portion  of  the  flow's  total  en¬ 
thalpy.  Chapter  7  details  these  effects  more  thoroughly.  Finally,  region  C  contains  cells  which 
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adapt  due  to  contributions  from  all  indices.  Note  that  an  actual  map  (e.g.  Ftg.  4.4)  contains  no 
cells  in  either  region  B  or  C. 


To  address  the  question  of  which  cells  adapt  due  to  which  parameters,  consider  the  Fig 
4.6  sequence  below. 


FIGURE  4.6 

Adaptation  baaed  on  different  parameter! 

These  figures  contain  grids  for  hot  nitrogen  flowing  at  Mach  5.66  over  a  blunt  wedge. 
Since  the  freestream  is  2390  K,  considerable  chemical  activity  occurs.  Using  p  as  an  adaptive 
index  captures  the  expansion  fan  and  developing  entropy  layer,  while  concentration  differences 
identify  the  chemical  relaxation.  Due  to  the  strong  coupling  as  the  flow  tends  toward 
equilibrium  at  the  stagnation  point,  both  indices  define  this  region. 

Threshold  Selection 

Often  great  debate  surrounds  the  topic  of  what  constitutes  a  "feature."  This  may  result  in 
elaborate  schemes  for  threshold  selection  (see  for  example  (7)  and  (33)).  In  practice,  however, 


the  question  is  not  critical  and  a  relatively  broad  range  will  adequately  capture  flow  physics. 
Different  thresholds  adapt  different  numbers  of  cells,  but  both  can  capture  the  feature 
reasonably  well  and  small  changes  in  "good"  threshold  value  should  not  radically  alter  the 
adaptation  pattern. 

As  an  example,  increasing  or  decreasing  Tadapt  by  10%  from  the  0.75  level  in  Figure  4.4 
would  have  little  effect  on  the  ultimate  number  of  cells  divided.  However,  if  located  within  the 
"crotch"  in  the  map,  slight  changes  in  value  would  produce  large  changes  in  the  size  of  adapted 
regions.  In  fact,  were  Tadapt  reduced  to  below  about  0.4,  the  arc  would  cut  arbitrarily  through 
the  data  based  on  meaninglessly  small  differences  in  adaptive  indices.  The  resulting  embedded 
grid  would  suffer  not  only  from  ragged  interface  patterns  but  also  from  problems  with  islands 
and  voids.  Threshold  values  for  shock  fit  problems  typically  range  from  0.6  -  0.7  since  the 
shock  fitting  has  eliminated  freestream  cells,  while  those  for  approaches  that  include  ffeestream 
cells  are  smaller  by  approximately  a  factor  of  two. 

Figure  4.7  traces  an  adaptive  sequence  cm  a  Mach  12  flow  over  a  blunt  wedge  in  air  at  an 
altitude  of  45  km.  Oxygen  dissociates  noticeably  but  the  stagnation  enthalpy  is  insufficient 
to  dissociate  appreciable  amounts  of  nitrogen.  These  maps  plot  normalized  differences  in 
oxygen  mass  fraction  against  normalized  density  differences.  For  this  2-D  example,  the  body 
size  is  such  that  the  reaction  completely  relaxes  in  about  1%  of  the  nose  radius,  producing  very 
stiff,  near  equilibrium,  behavior  near  the  shock. 

The  first  part  of  the  sequence  (4.7a)  shows  the  adaptation  map  for  the  original  grid  on  the 
left  and  the  adaptation  resulting  from  Tadapt  =  0.65  on  the  right  Since  the  chemical  behavior  is 
severe,  only  1  row  of  cells  extends  along  the  vertical  axis,  these  cells  corresponding  to  the 
string  of  data  extending  vertically  on  the  map.  Additionally,  the  scheme  adapts  cells  containing 
large  density  differences.  By  removing  the  extremum,  this  process  effectively  spreads  out  the 
remaining  pile  of  cells.  After  re-converging  the  solution  two  orders  of  magnitude,  the  adapta- 


tion  map  of  4.7b  is  essentially  a  close-up  of  the  original  map  and  reexamines  the  field  after  one 
level  of  embedding.  Notice  that  the  shock  cells  still  span  the  largest  chemical  gradients  -  i.e. 


Figure  4.7 

Adaptive  sequence  on  a  blunt  note. 

the  reaction  is  still  buried  within  this  first  set  of  cells.  The  differences  in  density,  however, 
exhibit  more  scatter,  are  more  evenly  distributed  across  the  abscissa,  and  a  few  cells  still 
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contain  large  differences.  Subdividing  this  map  with  the  same  0.6S  threshold  value  results  in 
the  grid  to  the  right.  Here  we  see  two  levels  of  embedding  indicating  the  relaxation  zone, 
entropy  layer  and  expansion  fan. 

Figure  4.7c  depicts  the  map  after  converging  on  the  doubly  embedded  grid.  The  cells 
are  now  evenly  distributed  across  the  abscissa  and  none  still  contain  extreme  density 
differences,  demonstrating  that  the  density  varies  smoothly  in  computational  space.  Only  one 
line  of  cells  still  contains  severe  chemical  behavior,  indicating  that  the  majority  of  the  relaxation 
still  spans  only  one  cell.  Notice,  however,  that  this  string  of  cells  demonstrates  little  variation 
in  density.  The  adaptive  process  has  separated  the  adaptive  indices. 

The  last  two  plots  (Figure  4.7c)  show  the  adaptive  map  after  a  third  adaptation.  Here  the 
abscissa  is  completely  diffuse  and  only  shock  cells  contain  large  chemical  gradients. 
Examining  the  actual  chem'cal  behavior,  we  see  that  despite  these  large  differences  across  the 
first  cell,  the  grid  completely  captures  the  chemical  relaxation.  Moreover,  the  initial  linear 
decay  contains  little  if  any  third  order  truncation  error  making  further  adaption  pointless  (and 
expensive). 


4.3  Unstructured  Grids 

The  mechanics  of  adaptation  revolve  around  the  use  of  unstructured  grids.  The  addition 
of  adaptation  to  any  unstructured  solver  is  usually  straightforward.  An  unstructured  grid 
should  not  be  considered  a  2-D  net  of  cells,  or  even  as  nodes  with  typical  i,  j  addresses. 
Instead  the  grid  is  simply  a  collection  of  cells  with  some  explicitly  defined  "connectivity" 
relating  the  cells  to  the  physical  domain. 
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An  unstructured  solver  operates  on  each  cell 
independently.  This  fact  alone  rules  out  most  implicit 
integration  algorithms.  Ni's  (29)  scheme  works  well 
because  the  calculation  proceeds  on  a  cell  by  cell  basis  and 
still  provides  second  order  spatial  accuracy  at  the  nodes. 

Moreover,  the  line  integration  is  quite  tolerant  of  misshapen 
cells  and  retains  its  accuracy  for  reasonably  large  degrees  of 
skewing  and  stretching  (7). 

Applying  the  integration  scheme  of  the  last  chapter  to  the  cell  in  Figure  4.8  requires: 

/.  Geometrical  information  describing  size,  shape,  etc.. 

ii.  The  ability  to  locate  neighboring  cells  for  any  node 

iii.  Knowledge  of  nodes  which  define  the  cell's  boundaries. 

Each  node  receives  a  unique  node  number  in  lieu  of  a  more  typical  i,j  address  permitting  a  two 
element  geometrical  pointer  to  locate  each  node  in  physical  space. 

GEOM(  1 ,  node  number)  =  1st  coordinate  in  physical  space 
GEOM(2,  node  number)  =  2nd  coordinate  in  physical  space 


FIGURE  4.8 

General  cell  in  unstructured  domain. 


Since  the  calculation  proceeds  cell-by-cell,  we  need  to  determine  which  cells  surround  any 
node.  The  neighbor  pointer  contains  four  elements  for  each  node 

NEIB(1,  node  number)  =  SW  cell  number 
NEIB(2,  node  number)  =  SE  cell  number 
NEIB(3,  node  number)  =  NE  cell  number 
NEIB(4,  node  number)  =  NW  cell  number 

FIGURE  4.9 

Cell  numbers  surrounding  a  node. 


If  the  node  is  along  a  boundary,  certain  cells  may  not  exist  and  this  pointer  receives  a  zero 
value. 


Once  inside  any  cell,  the  cell  pointer  locates  the  nine  nodes  defining  its  boundaries  (Fig. 


4.8). 


ICELL(1,  node  number) 

= 

Central 

node  number 

ICELL(2,  node  number) 

= 

SW 

node  number 

ICELL(3,  node  number) 

= 

S 

node  number 

ICELL(4,  node  number) 

= 

SE 

node  number 

ICELL(5,  node  number) 

= 

E 

node  number 

ICELL(6,  node  number) 

= 

NE 

node  number 

:  ■'ELL(7,  node  number) 

ss 

N 

node  number 

ICELL(8,  node  number) 

= 

NW 

node  number 

ICELL(9,  node  number) 

= 

W 

node  number 

Again,  if  a  node  does  not  exist,  its  location  in  the  cell  pointer  receives  a  zero. 

These  computational  structures  make  clear  the  advantages  of  the  integration  scheme  as  pre¬ 
sented  in  the  previous  chapter.  By  formulating  "cell  changes"  and  "distribution  formulae,"  we 
may  apply  this  node-based  scheme  in  a  completely  unstructured  manner  on  general  moving  or 
stationary  grids. 

In  addition,  other  pointers  can  be  defined  for  convenience.  The  overall  system  was 
developed  by  Dannenhoffer  and  Baron  (6)  and  has  been  used  with  minor  modifications  by 
(17),  (33)  and  (34)  among  others. 


4.4  Cell  Division 

Several  authors  discuss  the  process  of  cell  division  (7,  33,  and  34)  in  some  detail.  In 
general,  these  algorithms  are  tailored  to  a  specific  data  structure  and  pointer  system.  After  first 
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describing  the  steps  which  all  such  algorithms  must  include,  this  section  discusses  the  present 
approach. 

Cell  division  not  only  adds  new  cells  and  nodes  to  the  data  structure,  but  also  changes 
the  global  grid  connectivity.  Obviously  the  pointer  system  must  account  for  these  changes. 
The  matrices  must  absorb  the  new  cells  and  nodes  seamlessly,  insuring  that  new  cells  require 
no  special  treatment  by  the  integrator. 

The  steps  which  all  cell  division  algorithms  must  include  are  ( Fig.  4.10): 

a.  Find  cell  number  of  the  tagged  cell 

b.  Find  cell  numbers  of  the  12  surrounding  cells  (zero  if  no  cell  exists) 

c.  Find  node  numbers  of  nodes  2-9  if  they  exist 

d.  Create  new  node  at  center  (Fig.  4. 10,  node  1) 

e.  Create  new  nodes  at  the  center  of  each  face  if  they  do  not  exist  already 

/.  Update  pointers  of  cell  being  divided  to  reflect  new  comer  nodes  and  new  size 

g .  Create  four  new  cells  by  appending  Cell  pointer 

h .  Update  neighbor  pointers  of  nodes  1-9 

i.  Inform  surrounding  cell  pointers  of  any  new  nodes 

j.  Find  boundary  location  for  new  boundary  nodes 

k.  For  moving  grids,  specify  position  of  new  nodes  relative  to  existing  nodes 
(see  Chap.  5) 

Since  the  state  vector  and  auxiliary  properties  (p,  T,  h0,  etc.)  reside  at  nodal  locations, 
new  nodes  must  be  accompanied  by  initial  guesses  for  these  variables.  Properties  at  new  face 
nodes  take  values  equal  to  the  average  of  the  comer  nodes,  while  properties  at  central  nodes 
take  the  simple  average  of  all  four  comer  nodes. 
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FIGURE  4.10 

Cell  A  tagged  for  division. 


Naturally,  new  nodes  which  lie  on  boundaries  require  special  attention.  In  addition  to 
finding  the  correct  position  of  the  new  node  (either  from  interrogation  of  the  surface  splines,  or 
from  metrics  of  surrounding  nodes)  the  boundary  pointer  must  receive  this  new  information. 
By  looping  through  the  new  nodes  after  division,  it  is  relatively  simple  to  update  the  pointer  to 
account  for  new  boundary  nodes. 

Implementation 

Implementing  these  steps  requires  specific  decisions  about  data  structure,  appending 
pointers,  preserving  parent  cell  information,  etc..  As  a  result,  implementation  differs  with  each 
investigator. 

By  avoiding  "holes"  in  the  data  structure  where  a  cell  "existed"  before  division  and 
placing  new  cells  and  nodes  at  the  end  of  their  respective  pointers,  the  procedure  described  here 
results  in  an  adapted  data  structure  indistinguishable  from  the  original  grid.  Thus,  the  solver 
never  need  know  when  it  operates  on  an  adapted  cell.  This  permits  sequential  integration  of 
each  cell  without  ever  asking  if  a  cell  has  been  divided  and  without  skipping  over  any  vacancies 
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in  the  pointer.  While  very  attractive  from  the  viewpoint  of  simplicity  and  elegance,  this  system 
does  not  retain  parent  information  and  does  exclude  the  possibility  of  applying  a  multigrid 
accelerator. 

For  an  arbitrary  tagged  cell  n  in  a  grid  of  M  cells  Figure  4. 1 1  details  the  cell  division 
process.  Cell  n  shrinks  in  size  requiring  new  connectivity  (ICELL,  and  NEIB  pointers).  The 
scheme  numbers  the  three  new  cells  counterclockwise  M  +  1  through  M  +  3.  These  new  cells 
are  "connected”  to  the  appropriate  nodes  (1  — >9)  by  appending  the  ICELL  pointer.  All  nodes 
then  reconnect  to  their  surrounding  cells  by  updating  the  NEIB  pointer,  completing  the  pro¬ 
cess. 


Before  Division  After  Division 

FIGURE  4.11 

Cell  division  and  pointer  updating. 

Note  the  simplicity  of  this  scheme.  The  updated  ICELL  pointer  has  three  additional 
entries  and  no  vacancies.  The  NEIB  and  GEOM  pointers  extend  by  the  number  of  new  nodes 
with  a  minimum  of  complexity. 

References  7,  33  and  34  offer  more  elaborate  treatments  of  the  cell  division  process. 
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4.5  Cell  Interfaces 


FIGURE  4.12 

Interface  in  computational  mesh. 

Grid  refinement  and  adaptation  introduce  internal  boundaries  between  cells  of  different 
size.  Figure  4.12  depicts  such  an  interface.  Here,  only  one  column  of  cells  actually  differs 
from  those  considered  in  the  previous  chapter.  The  integration  scheme  holds  equally  well  on 
the  fine  cells  at  the  right  and  the  coarse  cells  to  the  left  of  A  and  B.  Each  of  these  has  exactly 
four  nodes.  However,  cells  bordering  on  the  embedded  region  (A  and  B)  always  contain  at 
least  one  hanging  node  and  require  special  treatment.  To  gain  insight  into  the  scheme's 
behavior  near  such  boundaries,  examine  the  two  situations  depicted  below. 
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FIGURE  4.13 

Comparison  of  interface  cells  with  cells  of  a  uniform  grid. 

Consider  the  model  convective  equation 

U,  +  Fx  =  0 .  (4.4) 

with  the  cell  and  flux  changes  for  any  cell  defined  by 

and 

AF  =  Fu&U ' 

Here  the  subscripts  w  and  e  represent  average  quantities  on  the  west  and  east  faces  of  any  cell. 

Since  this  is  essentially  a  1-D  equation  discretized  on  a  2-D  grid,  properties  remain 
constant  along  vertical  lines.  Now  examine  the  difference  between  nodes  /  in  Parts  I  and  II  of 
Figure  4.13. 

Although  node  /  receives  four  changes  from  its  surrounding  cells  in  both  cases,  the  two 
pictures  are  not  equivalent  Since  Ni's  scheme  is  simply  central  differencing  (on  cartesian 
grids),  we  expect  a  stretching  error  in  the  mesh  of  Case  II.  To  demonstrate  this,  examine  the 
changes  to  node  /  for  Case  I  and  Case  n,  5(7/  )i  and  8£//  )u  . 
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Case  I: 


By  Equation  (3.9): 


A  Uc 
4 


(4.5) 


Recalling  the  governing  Equation  (4.4),  and  considering  i-D  convection  only,  notice  that 
since  all  cells  are  identical,  8 UiE  =  5 U\F  =  8(7^=  8(/,f  .  Also,  the  lack  of  variation  in  the 

vertical  direction  makes  the  changes  at  node  i  from  B  and  C  equivalent  to  those  from  C  and  D  at 

/. 


For  the  first  case,  the  total  change  to  node  /  is  the  sum  of  the  contributions  from  the  four 
surrounding  cells.  Assuming  a  constant  CFL  number  defined  by 

A tdF 


CFL  s 


AxdU 


gives  Equation  (4.6). 


S{/,\  =  Suijc^F = 1{( AUc  +  AUdII  -  CFL)  +  (AUF  +  AUeH  *  CFL)] 


Case  II: 


Again,  using  Equation  (3.9)  gives: 


8  = 

BU„  =  Wd 

'-^>1 

8 

8  U,h.^L 

(4.7) 


where  WiG  =  WiH  and  8(//c  =  hUiD. 
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Summing  gives  the  total  change  at  node  /  across  the  interface. 

8(//)n  =  5 U i )c+d+g+h  =  c  +  AC/^Xl  -CFL)  +  (A Uc  +  AC7//X1  +  CFL)] 


We  may  now  compare  the  changes  at  node  /  between  the  two  cases.  Since  Cells  C  and  D  are 
identical  in  both  cases,  Equations  (4.6)  and  (4.S)  combine  into 

St//) n  =  &U,\  +  +  A  UE)  +  (A  UG  +  A  UH)]  (4  9 

The  term  in  brackets  represents  the  error  resulting  from  the  stretched  difference  stencil  in  the 
second  case.  Notice  that  this  first  order  error  disappears  when  A Uf  =  A Uq  and  A Ue  -  Af/// ; 
i.e.  when  the  difference  stencil  is  not  stretched. 

Having  examined  the  expected  error  due  to  grid  stretching  across  interfaces,  we  may  now 
take  up  the  issue  of  hanging  nodes  on  these  boundaries.  Referring  again  to  Figure  (4.13),  no¬ 
tice  that  cell  G  contains  a  hanging  node  /  on  its  eastern  face.  As  an  interface  node,  i  must 
suffer  from  the  same  stretching  error  as  /.  However,  we  hope  to  avoid  any  additional  error 
caused  by  distributing  incorrectly  from  cell  G  to  node  i.  The  scheme  should  treat  nodes  i  and  / 
identically. 


With  these  thoughts  in  mind,  the  net  change  to  node  i  in  Case  II  becomes: 

8f//)n  =  8(7  8  £/»)a  +  5(4:  +  Mt/U+  */8C//)c 

=  |(A  UB  +  A£/CX1  -  CFL)  +  m  +  CFL) 


(4.10) 


ku  and  ki  weight  the  contributions  to  the  upper  and  lower  nodes  on  the  eastern  face  in  order  to 
apply  them  directly  to  i.  For  example  if  ku  =  ki  -  V2  Equation  (4.10)  would  distribute  the  av¬ 
erage  of  Wuq  and  W[G  to  /. 


In  the  first  case 
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(4.11) 


5  Ui\  =  SUi)B  +  St/.lc  +  5t/i)tt  +  bUi)F 

=  t{{AUB  +  At/cXl  "  CFL)  +  (A^  +  AUpfl  +  CFL)] 

and  substituting  this  expression  into  (4.10)  for  St//)i  gives  (4.12). 

St4i  -  +  +  k,)-(AUA  +  4Uf)] 


However,  since  properties  vary  in  x  only,  A Ua  =  A UF  =  A Ue  and  A Ug  =  At///,  giving: 


8£/,)n  =  St/A  +  ii±fS=!{-(At/F  +  AUe)  +  A£/0(*.  +  *,)]  (4  n) 

Treating  nodes  i  and  /  identically  requires  that  the  change  at  i  equal  the  change  at  /.  Comparing 
(4.13)  with  that  for  node  lu  and  recalling  that  8(7, )i  =  8t//)i  leads  to: 

AUG(ku  +  ki)  =  (At/c  +  At/*) 


But,  since  A  Ug  =  At///, 

A  Uciku  +  ki)  =  2At/c 
or  ku  +  ki  *  2 

Since  no  reason  exists  favoring  either  node  u  or  / 

ku  =  k[  =1 


(4.14) 

(4.15) 


This  result  implies  that  the  mid-face  node  i 
should  receive  the  sum  of  the  changes  dis¬ 
tributed  to  the  vertices  of  the  face.  Figure 
4.14  displays  this  interface  treatment  pictori- 
ally. 


It  is  worth  noting  that  this  result  contrasts  that  presented  by  Dannenhoffer  (7)  where  he 
suggested  ku  =  ki  =1/2-  More  recently,  Kallinderis  (18)  examined  several  interface  treatments 
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showing  Equation  (4.15)  to  be  not  only  conservative,  but  also  a  time  accurate  interface 
treatment. 

Actually,  the  case  of  ku  =  ki  =!/2  bears  special  attention.  If  F  varies  linearly  and  Fy  is 
constant,  the  cell  change  in  a  large  cell  like  G  would  be  exactly  twice  that  of  a  small  cell  (E,  F, 
etc.).  Re-examining  Eq.(4.13)  now  reveals: 

6  £/,) d  =  S  J,  +  +  A U^k,  +  *,)  ]  (4  , 6) 

but  since  A  Uia,%e  is  twice  AU small,  the  bracket  vanishes.  ku  =  k[  =V2  exactly  corrects  for  the 
grid  stretching  error  in  this  special  case.  By  re-weighting  the  difference  stencil,  it  removes  the 
inaccuracy  at  node  i.  Unfortunately,  it  is  not  so  easy  to  counteract  the  stretching  effects  at  u  or 
/  resulting  in  an  non-uniform  treatment  of  the  interface  before  reaching  the  steady  state.  Notice 
also,  that  this  fortuitous  canceling  within  the  bracket  of  Eq.  (4.16)  occurs  only  under  the 
special  conditions  of  Fy=  Const  and  linear  variation  of  F. 

Conservation  vs.  Accuracy 

As  noted  earlier,  the  interface  treatment  described  in  the  preceding  section  is  conservative, 
but  the  stretched  difference  stencil  leads  to  a  first  order  error  as  noted  after  Eq.  (4.9).  It  is 
reasonable  to  look  for  an  interface  treatment  which  is  both  conservative  and  accurate. 
Kallinderis  examines  this  question  at  length  and  discovers  that  the  second  difference  smoothing 
induces  a  first  order  mass  flow  error  in  accurate  treatments  (18,  pp.145).  Since  the  shock  layer 
may  be  thought  of  as  a  control  volume,  any  induced  mass  flow  from  non-conservation  will 
create  an  error  in  shock  stand-off  distance.  Thus  we  opt  for  conservation  over  accuracy  to 
preserve  shock  location  and  shape. 
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4.6  Integration  Scheme  on  Interface  Cells 


The  cell  in  Figure  4.14  contains  a  hanging  node  on  its  eastern  face.  Failure  to  include 
information  from  this  node  when  computing  the  flux  integral  around  such  a  cell  will  obviously 
result  in  a  conservation  error  across  that  face.  Such  an  error  is  unacceptable,  as  emphasized  in 
the  last  section. 

The  integration  scheme  of  Chapter  3  (or  those  in  Appendix  A)  accommodates  this 
hanging  node.  With  the  distribution  algorithm  already  discussed,  only  the  calculation  of  cell 
changes  requires  attention. 


Consider  the  more  general  cell  of  Figure  4.15.  When  computing  the  flux  integral 
Equation  (3.13)  the  flux  through  the  east  face  may  be  split  into  two  parts. 

» 

>  F dy  -  Gdx 

(M±I^ne  -  yE)  .  -  X£j 

(4.17) 

FIGURE  4.15 
General  interface  cell. 


Here  the  first  bracket  tracks  the  flux  through  face  1  while  the  second  follows  flux  through  face 
2.  Since  all  divided  cells  have  straight  faces,  massaging  Equation  (4.17)  results  in  the 
parabolic  form  below. 


ft  — 


.  |Gw*2CttG«)(v>g.<gj 


(4.18) 
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Implementation 

If  the  mid-face  node  is  equal  to  the  average  of  its  comer  nodes,  Equation  (4.18)  reduces 
exactly  to  Eq.  (3.14).  This  property  suggests  a  general  interface  treatment  for  interface  cells. 
When  the  solver  first  examines  an  interface  cell,  it  creates  virtual  nodes  at  the  midpoint  of  each 
face  not  already  containing  a  real  node.  By  averaging  comer  values  of  F,  G  and  W  at  these 
virtual  nodes,  cells  with  hanging  nodes  on  any  (or  all)  faces  may  be  treated  by  a  single 
integration  formula  -  one  that  expects  nodes  on  all  faces.  This  approach  avoids  the  expensive 
logic  statements  and  repetitive  subroutines  required  to  treat  all  possible  combinations  of 
hanging  nodes  separately.  Moreover,  interface  cells  may  be  treated  with  no  decision  statements 
whatsoever  by  defining  a  pointer  to  locate  which  faces  of  any  cell  contain  real  midface  nodes. 

In  2-D  cartesian  coordinates,  the  first-order  cell  change  for  a  cell  with  real  or  virtual 
nodes  on  all  faces  is: 

(F„„  +  2F„+Fre)  ►w  -  ym)  -  (C|T^*C|t|li.  -  wd 

(Bag*  + 5g)|y„-w>)-(fe  2<^w 

(F« -  -  Wwri  -  (G;g  *  -  xsw)  _  (4  19) 


AU)c  =  A/Wc  -y— i 
Ac 


Efficiency 

Despite  the  expense  of  defining  an  additional  pointer  and  creating  virtual  nodes,  the  result 
is  less  machine  code  than  the  decision  statements  and  subroutines  they  replace.  Also,  since 
nodal  placements  change  only  during  adaptation,  this  pointer  updates  only  during  adaptation 
steps  -  at  negligible  cost.  This  general  treatment  is  much  simpler  to  code,  and  since  interface 
cells  are  only  a  small  fraction  of  the  total  cell  count,  minor  savings  gained  by  more  elaborate 
treatments  appear  difficult  to  justify. 
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5 .  Unstructured  Shock  Fitting 


In  many  situations  it  is  reasonable  to  take  the  bow  shock  as  the  upstream  boundary. 
"Fitting"  the  shock  in  this  manner  forces  the  computational  domain  to  be  coincident  witn  the 
disturbed  physical  domain.  The  shock's  position  is  unknown  a  priori,  and  its  final  position 
depends  upon  the  physics  within  the  shock  layer  and  the  geometry.  This  underscores  the 
importance  of  formulating  the  governing  equations  in  strong  conservation  law  form,  and 
treating  interfaces  conservatively. 

In  practice,  most  blunt  body  flows  lend  themselves  to  bow  shock  fitting.  Nevertheless, 
at  some  altitude  the  Knudsen  number  becomes  sufficiently  large  to  prohibit  accurate  modeling 
with  a  discontinuous  treatment  of  the  shock  boundary.  Nominal  shock  thickness  estimates 
are  between  five  and  seven  mean  free  paths  thick.  For  example,  at  an  altitude  of  78  km  in  the 
standard  atmosphere,  the  Knudsen  number  based  on  freestream  conditions  and  a  one  meter 
nose  is  0.0028.  The  shock  is  approximately  0.01  R, ,  thick.  Taking  this  as  a  rough  upper  limit 
for  the  assumption  of  a  discontinuous  shock,  we  wish  to  examine  other  assumptions  behind 
the  current  modeling.  The  sustained  flight  corridor  of  Figure  1.1  shows  that  airbreathing  flight 
at  upper  altitudes  requires  Mach  numbers  large  enough  to  invalidate  the  relatively  simple 
chemical  model  described  in  the  second  chapter.  The  assumptions  within  the  gas  model  break 
down  before  the  shock  thickness  is  appreciable.  Thus,  for  the  class  of  sustained  flight 
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problems  correctly  treated  by  the  gas  model,  the  bow  shocks  are  thin  enough  to  warrant  shock 
fitting. 

Many  early  shock  fitting  approaches  relied  on  a  time  dependent  technique  developed  by 
Moretti  (27).  This  approach  casts  the  governing  equations  in  characteristic  variables  to  solve 
the  compatibility  equation  along  characteristic  lines.  Matching  the  results  with  those  from  the 
moving  shock  Rankine  -  Hugoniot  relations  determines  a  unique  shock  speed  and  results  in  a 
time-accurate  procedure. 

Such  an  approach  suffers  from  two  drawbacks  in  the  context  of  real  gas  problems  solved 
on  unstructured  finite  volume  meshes.  The  Reimann  invariants  do  not  remain  constant  through 
the  entropy  gradients  downstream  of  a  curved  shock.  In  perfect  gases,  this  problem  alone  is 
not  limiting  since  one  may  always  place  nodes  close  enough  to  the  shock  that  it  "appears"  lo¬ 
cally  planar.  However,  in  chemically  excited  flows,  the  rapidly  relaxing  region  just  behind  the 
bow  shock  is  strongly  non-isentropic  along  streamlines.  The  closer  the  cells  lie  to  the  shock, 
the  steeper  the  chemical  gradients  and  the  more  severe  this  situation  becomes.  Secondly,  it  is 
not  clear  that  explicit  forms  of  the  characteristic  variables  and  compatibility  equation  exist  for 
general  equations  of  state. 

Blottner  and  Larson  (3)  recently  proposed  a  shock  fitting  technique  that  avoids  those 
problems.  The  next  section  discusses  their  basic  technique  applied  to  ncnequilibrium  flow. 
While  much  simpler  in  both  concept  and  implementation,  the  technique  is  not  time  accurate. 

Despite  the  fact  that  the  flows  under  consideration  are  out  of  equilibrium,  the  frozen 
shock  Rankine  -  Hugoniot  relations  still  apply.  The  translational  mode  equilibrates  much 
faster  than  chemical  modes.  Typically  the  shock  builds  in  fewer  than  ten  collisions,  while 
thousands  occur  before  the  chemistry  reaches  equilibrium.  Thus,  for  the  present  gas  model,  it 
is  reasonable  to  treat  the  classical  shock  as  a  sharp,  frozen  discontinuity,  while  still  permitting 
nonequilibrium  chemistry. 


5.1  Perfect  Gas  Shock  Fitting 


After  first  integrating  the  entire  domain,  shock  boundary  nodes  receive  suggested 
changes  from  their  neighboring  cells.  The  shock  fitting  scheme  corrects  the  state  vector  at  such 
points  consistent  with  the  appropriate  moving  shock  Rankine  -  Hugoniot  relations  and  moves 
the  shock  in  response  to  the  developing  interior  solution.  The  interior  nodes  then  undergo  a  re¬ 
mapping  between  the  upgraded  shock  location  and  body  surface. 

Beginning  with  the  Rankine  -  Hugoniot  shock  jump  relations,  this  sections  develops  the 
shock  fitting  technique  with  an  eye  toward  implementation.  For  some  known  ffeestream 
normal  Mach  number  component,  the  normal  shock  equations  determine  the  state  vector  U 
downstream  of  the  discontinuity.  However,  if  the  shock  is  moving,  one  must  first  find  its 
speed  before  computing  U  (Figure  5.1).  In  a  reference  frame  moving  with  the  shock,  the 
shock  "sees"  a  relative  frees tream  Mach  number  whose  normal  component  is  Mrk  which 
results  from  the  vector  subtraction  of  the  absolute  normal  Mach  number,  Mnoo,  and  the  non- 
dimensional  shock  speed,  bn.  Mnoa  is  defined  as  the  normal  ffeestream  Mach  number  in  the 
body's  frame  of  reference  -  with  no  relative  shock  motion.  Figures  5.1  and  5.2  define  these 
quantities  pictorially. 

Mr„  =Mnm-bn  / c  |  \ 
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FIGURE  5.1 

Definition  of  normal  relative  Mack  number. 


FIGURE  5.2 

Definition  of  absolute  normal  Mach  number. 


To  determine  Mrk  we  employ  the  solution  from  the  most  recent  integration  step.  Pre¬ 
sumably,  this  step  produced  changes  at  the  shock  nodes  resulting  from  the  arrival  of  waves 
which  propagate  throughout  the  domain.  Over  the  course  of  convergence  this  information  ul¬ 
timately  determines  the  shock  stand-off  distance. 

The  same  information  provides  an  indication  of  shock  speed.  For  example,  at  the  shock 
node  1,  the  ratio  of  the  downstream  temperature  T\  to  the  freestream  temperature  Too  could  be 
interpreted  as  resulting  from  some  shock.  Such  a  temperature  ratio  indicates  a  "shock"  of 
specific  strength,  with  a  specific  associated  Mach  number.  In  this  way,  ri/r«o  defines  a 
specific  relative  Mach  number.  However,  the  basis  of  shock  strength  need  not  be  temperature 
ratio.  Ratios  of  p,  e,  p,  etc.  also  determine  equally  valid  relative  Mach  numbers.  It  is  this 

non-uniqueness  that  makes  time-accurate  re-formulation  of  this  procedure  unclear. 

« 

Since  density  is  a  state  variable,  the  density  ratio,  Pi/poo,  was  used  to  determine  in 
the  present  work.  However,  exploratory  investigations  using  pressure  or  temperature  ratios 
displayed  negligible  differences  in  convergence  behavior. 

The  value  of  Mrh  then  determines  all  properties  "just  behind"  the  shock  through  the 
shock  jump  relations  and  equation  of  state.  The  updated  solution  vector  applies  to  a  moving 
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shock  and  to  the  steady  state  as  the  shock  speed  vanishes.  Rewriting  (5.1)  for  the  shock 
speed: 


bn  =Mn-~MR. 


(5.2) 


Complete  Equations 

Figure  5.2  describes  the  situation  in  a  2-0  or  axisymmetric  domain.  Here,  note  that  the 
outward  normal,  n,  does  not  necessarily  correspond  with  the  £  direction,  indicating  that  this 
formulation  includes  non-orthogonal  coordinates.  The  density  ratio  after  the  integration  step 
determines  Afy?„  from  the  shock  jump  relations. 

Ml.- 2 


pi 


15.3^ 


If  the  shock  were  stationary  with  respect  to  the  body,  the  absolute  normal  Mach  number  Af„w 
would  be:  (in  a  reference  frame  attached  to  the  body) 


Mn-  =  A*-— 


(5.4) 


The  assumed  density  ratio  and  relative  Mach  number  now  completely  specify  conditions  behind 
the  shock.  The  x  and  y  velocity  components  are: 


v,-^Lf,+v. 


(5.5) 


P-MR'gtf ' 


and  the  local  speed  is  V2  =  u2  +  v2  .  Here  gtf  corrects  for  any  non-orthogonality  of  the  local 
&  C  transformation.  Rom  vector  calculus 


(5.6) 
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Since  the  shock  is  assumed  frozen,  the  equation  of  state  for  a  mixture  of  perfect  gases  deter¬ 
mines  the  temperature  downstream  of  the  shock. 


Pi 


(5.7) 


The  energy  expression  from  (Eq.  2.8)  completes  the  state  vector  behind  the  shock. 

Note  that  the  state  vector  along  the  shock  is  consistent  with  a  moving  shock,  and  ho  will 
vary  until  the  shock  stops  moving.  At  that  time,  ho  takes  on  its  freestream  value. 

In  a  time  At,  the  shock  moves  a  distance  Ac,  Ay  in  accordance  with  its  normal  shock 
speed  from  Equation  (S.2). 

(5g) 

Ayt  =  -AtbJjLgft 

In  these  equations,  the  node  time  step  is  simply  an  average  of  the  time  steps  from  the  two 
neighboring  cells  (e.g.  in  Figure  5.2  A/,-  =  Ar*  +  A ffl). 


5.2  Equilibrium  Shock  Fitting 

As  an  upper  limit  for  finite  rate  chemistry,  equilibrium  provides  a  convenient  basis  for 
comparison  with  very  rapid  chemistry.  Given  this  motivation,  we  now  develop  a  separate 
method  for  fitting  equilibrium  shocks.  The  discussion  revolves  around  a  specific  equilibrium 
gas  model,  but  applies  to  any  gas  in  chemical  equilibrium. 
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For  simplicity,  consider  the  reaction  A2  +  M  2A  +  M  as  described  by  Lighthill's  ideal 
dissociating  gas  model.  In  equilibrium,  the  shock  jumps  are  accompanied  by  changes  in  y 


across  the  shock.  For  Lighthill's  model, 

y  —  4  ±  Ql 

1  3  (5.9) 

Here  a  is  the  degree  of  dissociation  (36,  pp.  159-161).  Thus  the  ratio  of  specific  heats 
depends  upon  the  degree  of  dissociation.  However,  simply  changing  y  in  the  perfect  gas 
Rankine  -  Hugoniot  relations  neglects  the  aSd  energy  absorbed  in  dissociation  (see  for 
example  Eq.  2.6).  Obviously,  the  resulting  shock  would  not  preserve  stagnation  enthalpy,  and 
would  violate  energy  conservation. 

Ruling  out  the  use  of  closed  form  shock  jump  relations  leaves  the  need  to  perform 
equilibrium  normal  shock  calculations  and  then  find  the  correct  shock  speed. 

Mechanics 

Referring  again  to  Figure  5.1,  the  flowfield  solution  from  the  previous  iteration  provides 
a  pressure  jump  across  the  shock.  Pi /poo.  Using  the  equation  of  state  with  an  assumed  value  of 
a\  gives  the  density  just  downstream  of  the  shock. 

_  pi _ 

Pl  RAtTil+ai)  (5.10) 

Conservation  of  momentum  through  the  discontinuity  results  in  an  expression  for  the  relative 
normal  velocity  at  1  (36,  pp.  179). 

fjl  _  Pi  (Pi 

Rm  P-(Pi-P~>  (5.11) 


Energy  conservation  determines  a  corresponding  value  of  enthalpy  behind  the  shock. 


h\  -  /u  + 


(5.12) 


88 


Recalling  the  enthalpy  expression  for  the  ideal  dissociating  gas: 


h=RA2[{4  +  a)T  +  a  ed] 


Re-arranging,  expresses  the  temperature  downstream  of  the  shock  as  a  function  of  enthalpy 
and  mass  fraction. 


Ti  = 


1 

’  h\ 

(4  +  ax) 

Ra2 

-a\Sd 


(5.13) 


and  this  temperature  permits  a  better  estimate  of  the  atomic  mass  fraction  at  1. 

a  _  ,L±h\+  4 k  where  k  a  — e®^ 

2  P  (5.14) 

However,  (5.13)  assumed  the  previous  value  of  a,  and  we  must  now  compute  a  new 
temperature  consistent  with  the  mass  fraction  determined  by  (5.14).  Iteration  between  Equa¬ 
tions  (5.13)  and  (5.14)  determines  a  unique  mass  fraction  and  temperature.  Then  (5.10)  pro¬ 
vides  a  consistent  value  of  the  density  ratio  across  the  shock  and  (5.1 1)  yields  an  improved 
shock  speed. 


This  doubly  iterative  process  initially  may  appear  computationally  expensive.  However, 
it  is  applied  only  to  shock  points,  and  the  net  cost  remains  small.  After  obtaining  a  converged 
relative  Mach  number,  the  shock  speed  and  subsequent  motion  follow  as  in  the  previous 
section. 


Some  Notes  on  Equilibrium  Shock  Fitting 

The  method  described  above  is  not  elegant  One  possible  alternative  might  make  use  of 
the  fact  that  thermodynamic  equilibrium  minimizes  the  Gibbs  free  energy.  Setting  up  the 
minimization  problem  permits  the  use  of  LaGrange  multipliers  constrained  by  the  existing 
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pressure  or  density  jump.  Certainly  extending  the  model  to  include  multiple  coupled  reactions 
(NO  present)  with  many  reaction  paths  would  demand  a  more  elaborate  approach. 

As  a  final  comment,  the  stiff  transcendental  function  in  Equation  (S.14)  results  in  extreme 
sensitivity  of  a  to  small  changes  in  T.  Specifically,  a  tends  to  oscillate  between  0.0  and  1.0, 
hampering  convergence  of  this  inner  iteration.  A  simple  fix  is  to  reduce  changes  in  a  by 
approximately  70%,  thus  decreasing  changes  to  T,  and  resulting  in  a  more  tractable 
convergence  behavior. 


5.3  Moving  Unstructured  Meshes 

A  moving  inflow  boundary  implies  that  the  grid  must  redistribute  itself  to  properly  track 
the  shock's  motion.  For  structured  grid  calculations,  this  merits  relatively  little  attention. 
However,  in  the  realm  of  the  current  unstructured  domain,  more  care  is  required. 

The  grid  redistribution  requires  knowledge  of  shock  position  and  information  about  a 
particular  node's  history  in  the  grid  (stretching  or  skewing).  Since  all  previous  discussions 
involve  current  cell  information  only,  such  global  location  requests  do  not  fit  within  the 
framework  of  an  unstructured  scheme. 


Reference  points  for  moving  nodes  on  the  base  (original)  grid. 


The  endpoint  pointer  fixes  all  nodal  locations  with  respect  to  any  set  of  four  reference 
nodes.  The  pointer's  six  elements  store  reference  node  numbers  as  well  as  x  and  y  position 
with  respect  to  these  references.  Referring  to  the  nodes  on  the  unadapted  mesh  of  Figure  5.3, 
the  pointer  fills  as  follows. 

ENDPT(1,c)  =  body  node  a) 

ENDPT(  2,  c )  =  body  node  Hb) 

ENDPT(  3,  c )  =  %  x  of  total  distance  IF  measured  from  a 
ENDPT(  4,  c )  =  %  y  of  total  distance  aF  measured  from  a 
ENDPT(  5,  c )  =  body  node  tt  (a) 

ENDPT(  6,  c )  =  body  node  U  (b) 


These  definitions  locate  original  grid  nodes  by  their  position  along  lines  of  constant  £.  Since  it 
is  based  upon  the  original  grid,  this  pointer  preserves  any  curvature  or  clustering  of  grid  lines. 


After  each  adaptation,  every  new  cell  must  have  faces  matching  those  of  its  neighbors 
across  cell  interfaces.  Failure  to  keep  these  faces  linear  results  in  higher  order  conservation 
errors  across  the  interface  when  integrating  adjacent  cells.  Unfortunately,  referring  mid-face 
nodes  to  shock  and  body  points  does  not  insure  linear  cell  faces  as  the  grid  moves,  and  we 
seek  an  alternative  solution. 
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Figure  5.4  examines  this  situation  in  more  detail.  If  midface  and  central  nodes  rely  on 
body  and  shock  locations,  there  is  no  guarantee  that  linear  interfaces  will  be  preserved.  To 
avoid  this  problem,  the  endpoint  pointer  positions  new  nodes  by  the  comers  of  their  supercells. 
For  the  example  shown  in  the  sketch,  all  new  nodes  should  rely  on  the  comer  nodes  of  the 
original  Cell  A  .  Specifically,  for  new  nodes  created  from  the  division  of  Cell  A,  the  ENDPT 
pointer  entries  1, 2, 5,  and  6  will  refer  to  comer  node  data  for  the  original  Cell  A.  Using  these 
reference  points,  each  midface  node  remains  half  way  between  the  comers,  falling  exactly  on 
the  interface.  Through  this  more  general  use  of  the  pointer,  the  re-mapping  algorithm  can  treat 
new  nodes  exactly  as  original  nodes,  efficiently  creating  a  new  grid  without  any  decision 
statements  (and  minimizing  computational  effort).  Finally,  note  that  this  pointer  permits 
remapping  grids  with  any  number  of  embedded  levels,  while  simultaneously  maintaining  all 
cell  interfaces  as  straight  lines. 


FIGURE  5.4 


II.  Incorrect  redistribution 


III.  Correct  redistribution 


Preservation  of  linear  interfaces  in  an  unstructured  domain  with  multiple  levels  of  adaptive  embedding. 


5.4  Behavior  of  Shock  Fit  Solutions 

As  an  initial  condition,  all  points  within  the  domain  (Fig  3.14,  including  shock  points) 
might  contain  ffeestream  values.  In  this  case,  no  shock  discontinuity  exists,  and  the  initial 
conditions  exactly  satisfy  the  governing  equations  at  all  points  not  on  the  body.  With 
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freestream  values,  all  shock  points  span  a  density  ratio  P»Vp u  *1.0  and  the  shock  degenerates 
to  an  acoustic  wave.  Confirming  this.  Equation  (5.3)  predicts  a  relative  normal  Mach  number 
of  unity,  and  (5.5)  through  (5.6)  predict  the  same  for  the  velocity  ratio,  temperature  ratio  and 
pressure  ratio  at  all  shock  points  and  the  "shock"  moves  inward  at  a  non-dimensional  speed  bn 
Of  Mm  — '  1 . 


0.0  1.0  20  ao  4.0  5.0 


Thousands  of  ttsrations 


FIGURE  5.5 

A  typical  history  of  stand-off  distance  with  iteration,  converged  from  a  uniform  freestream  initial 
condition.  (Mach  5  frozen  flow  over  a  axisymm eerie  nose.) 


During  the  first  time  step,  the  explicit  solver  transmits  knowledge  of  the  body  to  adjacent 
cells  one  removed  from  this  boundary.  Meanwhile,  the  shock  continues  progressing  towards 
the  body  at  a  speed  Af«  -  1.  On  a  grid  with  M  lines  between  the  shock  and  body,  the  shock 
will  not  respond  to  the  body's  presence  for  M  time  steps.  This  behavior  is  a  direct  result  of  the 
explicit  algorithm,  and  results  in  a  time  lag  between  the  shock  and  the  developing  solution. 
Moreover,  it  creates  a  feedback  loop  wherein  the  shock's  position  and  shape  determine 
properties  within  the  domain,  but  the  shock  cannot  respond  immediately  to  these  changing 
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properties.  In  this  respect,  the  problem  resembles  an  elliptic  calculation  (albeit  one  with  an 
internal  time  lag). 

Figure  5.5  shows  a  typical  history  of  shock  motion  recorded  while  converging  a  Mach  5 
solution  from  freestream  conditions.  Initially,  the  shock  resides  some  distance  from  the  body 
and  moves  toward  it  as  described  above.  Then,  as  pressure  waves  cross  the  layer,  the  Rankine 
-  Hugoniot  jumps  build,  slowing  the  shock  before  reversing  its  motion  altogether.  The  shock 
backs  away  from  the  body  and  eventually  comes  to  rest  at  some  final  stand-off  distance  deter¬ 
mined  by  the  interior  solution. 

Although  this  example  of  shock  adjustment  was  computed  from  a  special  (freestream) 
initial  condition,  the  shock's  behavior  results  more  generally  from  wave  propagation  in  the 
explicit  shock  fit  domain.  In  practice,  the  shock  adjusts  its  position  in  a  similar  manner  to  all 
disturbances  (changes  in  freestream  conditions,  chemistry,  etc.). 

Since  adaptation  selectively  adds  cells  to  the  domain,  it  changes  the  wave  propagation 
speed  through  various  regions  of  the  domain  (in  physical  space).  As  a  result,  adapting  during 
the  initial  gross  shock  motion  may  disrupt  the  shock's  shape  enough  to  destabilize  the  solution. 
Simply  delaying  until  the  stand-off  distance  stabilizes  avoids  this  shortcoming. 
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Ill  PRESENTATION  AND  DISCUSSION 
OF  RESULTS 


A  DEC  Microvax  3200,  nominally  rated  at  3  Million  Instructions  Per  Second  (MIPS),  provided 

processing  support  for  all  test  cases  in  this  work. 


6.  Physical  Phenomena  in  High  Temperature  Flows 
Over  Simple  Geometries 

The  gas  model  and  adaptive  integration  scheme  discussed  in  Parts  I  and  II  led  to  the 
development  of  a  series  of  adaptive  hypersonic  CFD  codes.  These  codes  permit  computation 
of  inviscid  and  viscous,  chemically  reacting  flow  over  2-D  and  axisymmetric  bodies.  The 
additional  development  of  a  2-D  identical  equilibrium  Euler  solver  provided  a  datum  for  flows 
demonstrating  very  small  departure  from  equilibrium. 

After  a  brief  validation  study,  focus  turns  to  the  question  of  using  adaptation  to  resolve 
physical  behavior  within  the  reacting  shock  layer.  Both  inviscid  and  viscous  results  support  a 
fundamental  discussion  of  chemical  length  scale  behavior  within  the  shock  layer.  Finally,  as 
we  consider  flows  with  higher  Mach  numbers,  interest  broadens  to  coupled  and  uncoupled 
multiple  reacting  mixtures,  examining  both  the  formation  and  impact  of  NO  within  a  gas  cap, 
and  the  degree  of  coupling  between  reactions.  These  insights  form  a  basis  for  then  evaluating 
the  effectiveness  of  adaptive  grid  embedding  in  hypersonic  shock  layers  in  Chapter  7. 

6.1  Basic  Examples  and  Algorithm  Verification 

Before  looking  at  the  detailed  physics  and  flow  phenomena  within  the  shock  layer, 
validation  studies  were  completed  to  lend  credibility  to  both  the  basic  algorithm  and  gas  model. 
This  section  outlines  three  such  studies,  designed  to  verify  different  aspects  of  the  solution 
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algorithm.  First  is  a  classic  perfect  gas  test  case,  which  demonstrates  both  the  accuracy  and 
conservation  of  the  Lax-Wendroff  algorithm  and  unstructured  shock  fitting  procedure.  A  more 
rigorous  test  e  amines  the  adaptation  methodology  and  gas  model  behavior  in  severe  reacting 
flows.  A  final  example  evaluates  the  effectiveness  of  the  present  method  in  tracking  multiple 
reaction  paths  and  coupled  reactions  in  a  multi-reaction  mixture. 

Perfect  Gas  Test  Cases 

Figure  6.1  displays  the  traditional  variation  of  stand-off  distance  with  Mach  number  for  a 
2-D  cylinder  exposed  to  a  crossflow  and  an  axisymmetric  hemispherical  cylinder  test  case. 
Reference  22  provided  the  comparison  curves,  and  the  data  show  numerical  experiments  at 
Mach  3, 5, 7.5  *nd  10  in  a  y=  1.4  perfect  gas. 


FIGURE  6.1 

Comparison  of  shock  Stand-off  with  Reference  22. 

These  data  appear  to  deviate  slightly  from  the  reference  curves  at  high  Mach  numbers. 
This  discrepancy  arises  as  a  direct  result  of  the  half  excited  vibrational  state  in  the  present  mod¬ 
eling.  The  normal  shock  relations  predict  a  density  jump  across  the  shock  of 
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(6.1) 


p.k_ 

P-  (y-l)Mi+2 

This  ratio  asymptotically  approaches  6  at  large  Mach  number  for  perfect  air  with  y  =  1 .4.  Since 
the  flow  behind  the  shock  is  subsonic  in  the  stagnation  region,  the  density  ratio  will  not  in¬ 
crease  much  as  the  flow  comes  to  rest  at  the  stagnation  point,  and  density  throughout  the  region 
remains  relatively  constant.  For  the  present  gas  model  with  its  half-excited  vibrational  state, 
kinetic  theory  predicts  a  frozen  ratio  of  specific  heats  of  73  in  the  absence  of  dissociation, 
which  leads  to  a  limiting  density  ratio  of  7.  With  approximately  the  same  mass  flux  entering  a 
shock  layer,  the  shock  in  4/3  test  cases  should  be  roughly  16%  closer  to  the  body  in  the  limit  of 
infinite  freestream  Mach  number.  At  finite  Mach  numbers,  however.  Equation  (6.1)  predicts 
slightly  less  pronounced  effects  due  to  y.  Gose  examination  of  the  figure  reveals  that  at  higher 
Mach  numbers,  the  calculated  stand-off  distance  approaches  this  proper  infinite  Mach  number 
limit.  For  example,  at  Mach  10,  the  axisymmetric  code  predicts  a  shock  location  13%  closer  to 
the  body  than  the  extrapolated  y=1.4  result 

A  second  classical  result  predicts  the  shape  and  position  of  the  sonic  line  within  the  shock 
layer  (13).  Figure  6.2  displays  the  sonic  lines  for  2-D  and  axisymmetric  calculations.  Both 
pictures  depict  Mach  5  flow  of  the  same  y=  V3  perfect  gas. 

The  sonic  line  in  the  2-D  flow  shown  on  the  left  clearly  demonstrates  the  high  curvature 
associated  with  sonic  line  behavior  above  Mach  2  in  2-D  supersonic  blunt  body  flows,  and 
displays  the  expected  acute  angle  with  the  body's  surface.  In  the  axisymmetric  case,  the  sonic 
line  takes  on  the  steeply  raked  profile  characteristic  of  3-D  blunt  body  flows.  The  third 
dimension  provides  an  additional  direction  for  the  flow  to  expand,  and  this,  combined  with  the 
thinner  shock  layer,  tends  to  flatten  the  profile.  As  is  typical  of  freestream  Mach  numbers 
greater  than  3,  the  angle  between  the  body  and  sonic  line  is  obtuse. 
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While  primarily  qualitative  evidence,  the  behavior  of  the  sonic  line  does  demonstrate  the 
Lax-Wendroff  scheme's  ability  to  model  the  physics  between  the  shock  and  body  •  even  in  the 
stagnation  region.  The  predicted  shock  stand-off  distances  lend  credence  to  both  the  shock  fit¬ 
ting  implementation  and  the  schemes  mass  conservation  properties. 


figure  6.2  Two-Dimensional  Axisymmenic 

Computed  Mach  contour*  of  Mich  5  frozen  flow  over  2-D  and  txisymmetric  circularly  blunted 
bodies  showing  location  of  of  the  sonic  line. 

Dissociating  Nitrogen 

While  effective  from  the  standpoint  of  understanding  and  completeness,  these 
comparisons  tend  to  be  primarily  qualitative.  Essentially,  the  above  results  demonstrate  that  the 
shock  fit,  Lax-Wendroff  technique  may  be  applied  to  supersonic  flows,  but  do  not  rigorously 
verify  the  technique's  accuracy,  nor  do  such  perfect  gas  computations  demonstrate  anything 
about  the  real  gas  model.  Therefore,  we  now  compare  the  current  method  with  both 
experimental  and  computational  results  published  by  M.  N.  Macrossan  (25). 

The  case  provides  a  very  stringent  test  of  the  scheme's  ability  to  model  high  temperature 
flows.  Table  6. 1  details  the  freestream  conditions  for  both  the  experiment  and  computation.  In 
both  situations,  hot  nitrogen  flows  at  approximately  Mach  5  past  a  15°  blunted  wedge.  The 
computational  conditions  were  chosen  to  match  momentum  flux  through  the  shock  and 
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stagnation  enthalpy  with  those  for  the  experiment.  Notice  that  test  conditions  are  "nominal” 
and  Macrossan  estimates  a  possible  2-7%  measurement  uncertainty  (25). 

Several  attributes  make  this  an  excellent  test  case.  First,  the  flow  over  much  of  the  blunt 
wedge  is  in  a  state  of  small  departure  from  equilibrium.  The  chemical  length  scale  behind  the 
shock  is  much  smaller  than  the  physical  dimensions  of  the  body.  Using  properties  behind  a 
frozen  normal  shock  for  the  freestream  conditions,  the  chemical  relaxation  length  is  slightly 
over  0.003  Rn •  In  terms  of  the  Damkohler  number  defined  at  the  end  of  the  first  chapter,  is 
greater  than  36,000.  Clearly,  we  expect  very  rapid  chemical  relaxation  near  the  shock.  More¬ 
over,  the  large  dissociation  energy  of  nitrogen  will  strongly  couple  these  chemical  effects  to 
other  flow  parameters  through  the  shock  layer.  This  very  strong  coupling,  combined  with  the 
very  rapid  dissociation  provides  a  severe  test  of  the  gas  model  since  over  40%  of  the  stagnation 
enthalpy  is  absorbed  by  dissociation  within  the  first  5%  of  the  shock  layer. 

TABLE  6.1 

Conditions  for  comparison  to  dissociating  nitrogen  test  case  (from  Reference  25). 


Real  shock  tunnel  Computation  with 

nominal  conditions  equilibrium  freestream 

Woo  (m/s) 

PooMm  3) 

Too  (K) 
o»(£3S.) 

Moo 

Rn  (mm) 

6.36  x  103 

4.41  x  10*2 

4415 

0.094  (frozen) 

4.55 

5.0 

6.31  x  103 
4.45  x  10*2 

5570 

0.064 

4.09 

4.91 

In  addition  to  the  stiffness  arising  from  the  disparate  time  scales  involved,  this  case 
demonstrates  the  behavior  of  the  gas  model  in  flows  with  small  regions  where  temperatures 
exceed  9,000  K.  Normal  shock  calculations  predict  temperature  ratios  of  about  4  across  the 
shock.  Since  the  freestream  temperature  is  approximately  5,000  K,  the  temperature 
immediately  behind  the  shock  is  in  excess  of  that  allowed  with  the  present  model.  However, 
the  rapid  chemical  relaxation  in  fact  quickly  absorbs  much  of  this  energy  and  drops  the 
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temperature  radically.  Most  of  the  flow  remains  within  the  range  of  applicability  of  the  half- 
excited  assumption. 

Finally,  capturing  the  relaxation  near  the  bow  shock  required  three  levels  of  embedding 
and  placed  interface  cells  in  the  stagnation  region.  This  demonstrates  the  ability  of  the  adapta¬ 
tion  to  resolve  flow  features,  while  preserving  conservation  across  interfaces.  Since  the  shock 
layer  may  be  thought  of  as  a  control  volume,  any  error  in  conservation  will  create  an  error  in 
stand-off  distance. 


Cunent  method  shown  in  lower  half. 

Figure  6.3  compares  density  contours  computed  with  the  present  model  with  those  from 
both  experiment  and  computations  performed  by  Macrossan.  Both  agree  quite  closely.  The 
figure  on  the  left  contains  experimental  results  above  and  results  from  the  adaptive  computation 
at  nominal  shock  tunnel  conditions  below.  The  comparison  on  the  right  displays  Macrossan's 
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solution  for  the  "equivalent  equilibrium  freestream"  above  a  result  from  the  current  method.  In 
Table  6.1  it  is  interesting  to  note  that  reference  (25)  intended  these  results  to  demonstrate  that 
the  "equivalent"  conditions  gave  identical  results  as  the  experiment,  and  considered  the 
flowfield  at  the  right  virtually  identical  to  that  at  the  left 


While  the  comparison  with  the  calculation  is  in  rather  good  agreement,  the  comparison 
with  experiment  differs  slightly.  The  density  in  the  flow  varies  from  roughly  10  to  2  and  the 
present  solution  shifts  the  2.5  and  3.0  contours  by  approximately  0.25pw.  This  points  to  a 


discrepancy  of  approximately  2.5%.  However,  its  worth  remembering  that  the  freestream 
conditions  held  an  uncertainty  of  approximately  5%,  and  therefore  the  comparison  is  within 


experimental  error. 


FIGURE  6.4 

Adaptive  computational  (rid  with  1640  nodes  for  distociating  nitrogen  flow  over  a  circularly 
blunted  15*  wedge  at  'equivalent'*  conditions  (Table  6.1). 


As  evidence  of  the  adaptive  scheme’s  ability  to  detect  and  resolve  flow  features,  consider 
Figure  6.4.  Here,  the  final  adapted  grid  for  the  "equivalent"  conditions  (Table  6.1)  case  clearly 
shows  the  relaxation  zone  and  expansion  fan  in  the  flowfield.  With  three  levels  of  adaptive  re¬ 
finement,  the  final  grid  contains  1640  nodes  or  slightly  more  than  half  the  3000  nodes  used  by 
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the  DSMC  calculation  of  the  reference.  The  RMS  momentum  residual  was  converged  two  or¬ 
ders  of  magnitude  between  adaptations.  The  adaptation  threshold,  7^^  was  set  at  0.35  based 
upon  examination  of  the  first  adaptation  map  and  held  at  this  level  for  the  remaining  two  adap¬ 
tations. 

The  source  did  not  estimate  the  Reynolds  number  or  other  viscous  parameters  in  die  hot 
nitrogen  test  gas.  A  rough  estimate  is  that  Reynolds  number  is  on  the  order  of  10s.  The 
interferograms  presented  in  Reference  25  seem  to  confirm  such  a  high  value  since  the  boundary 
layer  remains  quite  thin. 


We  now  focus  on  the  behavior  of  the  mixture  model  and  thermodynamics  as  presented 
in  Chapter  2.  Comparison  of  the  current  gas  model  is  made  with  two  other  models  described 
by  Candler  (5)  and  an  experiment  of  a  sphere  fired  in  air  at  5280m/s  performed  by  Lobb  (24). 
Candler’s  gas  models  include  both  a  single-  and  six-temperature  models.  The  latter  includes 
translational  temperature,  four  vibrational  temperatures  (W2. O2,  NO,  NO*),  and  an  electron 
temperature.  The  reference  has  demonstrated  that  the  six-temperature  model  predicts  the 
sphere's  stand-off  distance  and  shock  shape  to  within  experimental  accuracy,  and  therefore 
comparison  is  made  directly  to  the  six-temperature  model 


To  begin,  examine  the  stand-off  distance  implied  by  the  present  and  Ref.  5  density  ratio 
distributions  along  the  stagnation  streamline  in  Figure  6.5.  Table  6.2  details  assumed 
frees  tream  conditions. 
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TABLE  6.2 
Conditions  for 


In  figure  6.5,  the  gas  flows  from  left  to  right  and  the  stagnation  point  is  at  *1/5  =  0.0.  The 
present  computations  were  inviscid  and,  as  expected,  the  shock  stand-off  distance  differs  by 
about  5%  due  to  the  thickness  of  the  boundary  layer  (-10%)  at  this  low  Reynolds  number. 
Note  that  the  shock  fitting  in  the  present  method  results  in  a  discontinuous  shock  (as  would  be 
expected  at  these  conditions).  The  error  in  stand-off  distance  arises  directly  from  the  cool  ther¬ 
mal  boundary  layer  of  the  viscous  calculations.  There,  nearly  constant  pressure  and  lower 
temperature  force  density  to  increase  rapidly,  bringing  the  shock  closer  to  the  body.  In  the 
inviscid  portion  of  the  flowfield,  the  present  model  does  agree  within  plotting  accuracy  for  the 
more  elaborate  models. 


0.1  0.05  0 

T|^>,  Normalized  Distance  From  Stag.  Point 


FIGURE  6.5 

Stagnation  streamline  profiles  of  density  ratio  of  current  method  and  computations  of 
Candler  (5)  at  conditions  in  Table  6.2. 
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Figures  6.6a  and  6.6b  show  further  details  along  the  stagnation  streamline.  The 
temperature  profiles  display  especially  interesting  behavior.  Since  both  Ref.  5  models  use 
shock  capturing,  neither  resolves  the  temperature  peak  across  the  shock.  Those  results  for  the 
six -temperature  model  imply  that  vibration  remains  unexcited  across  the  translational  shock, 
supporting  the  hypothesis  of  a  frozen  shock,  and  that  significant  vibrational  nonequilibrium 
exists  throughout  much  of  the  shock  layer.  If  frozen,  the  temperature  peak  for  y=  1.4  should 
be  approximately  14,000  K.  The  present  y=  V3  model  has  a  peak  temperature  of  1 1,500  K 
just  downstream  of  the  shock.  The  one-temperature  model  (5)  allows  vibration  to  absorb 
energy  through  the  captured  shock,  and  results  in  a  10,000  K  peak  behind  the  shock.  As 
chemical  reactions  proceed,  the  curves  decay.  Both  single  temperature  models  decay  rather 
quickly  as  dissociating  O2  soaks  up  translational  energy.  The  multitemperature  model  decay  is 
less  rapid  since  its  dissociation  rate  is  linked  to  both  vibrational  and  translational  temperature. 


n/5,  Nocnalixad  Diajnca  From  SU|.  Poo*  llA  NonmBfrt  Diittne*  F«cm  S«»».  IVitat 

FIGURE  6.6a  FIGURE  6.6b 

Stagnation  streamline  profiles  of  ter-nerature  and  species  mass  fraction  of  current  method  and  computations  of  Candler 
(5)  at  conditions  in  Table  6.2. 

This  effect  also  explains  the  delay  (of  nearly  0.16)  experienced  before  the  onset  of  O2 
dissociation  for  the  six-temperature  model  (Fig.6.6b).  Although  the  "shock"  stands  at  0.1  Rn, 
O2  hesitates  until  0.09  Rn  before  dissociation  commences.  Candler  ascribes  this  behavior  to 
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nonequilibrium  within  the  vibrational  modes.  Otherwise,  the  figure  demonstrates  expected 
behavior  of  O2  and  NO  throughout  the  shock  layer.  Nevertheless,  the  precise  levels  of  NO 
disagree  amongst  the  models  and  the  inviscid  axisymmetric  code  used  in  this  comparison 
naturally  tracks  none  of  the  behavior  in  the  wall  boundary  layer. 

To  summarize,  although  the  current  modeling  differs  in  prediction  of  the  absolute  species 
concentration  levels,  it  demonstrates  the  technique's  ability  to  model  reaction  coupling,  with 
very  reasonable  prediction  of  the  gas  dynamic  variables. 

Note  that  the  temperature  scale  in  Figure  6.6a  reads  in  Kelvin.  The  majority  of  this  shock 
layer  is  above,  or  near,  the  9,000  K  limit  indicated  by  the  derivation  in  Chapter  2.  Despite 
this,  the  predicted  level  of  NO  molecules  remains  comparable  to  those  from  the  other  models. 

6.2  Effects  of  Length  Scale  Behavior 

The  temperature  relaxation  seen  in  Figure  6.6  characterizes  nonequilibrium  shock  layer 
calculations.  The  ratio  of  the  chemical  relaxation  length  to  a  relevant  physical  reference  length, 
XlLref'  is  a  most  important  parameter  in  the  description  of  nonequilibrium  situations.  Such 
parameters  characterize  the  relative  importance  of  the  chemical  modeling  and  determine  the 
distance  required  for  the  temperature  to  decay,  or,  equivalently,  for  the  density  to  increase.  As 
a  direct  result,  the  shock  layer  thickness  depends  strongly  upon  the  degree  of  nonequilibrium. 

General  Length  and  Time  Scale  Behavior  in  Inviscid  Hypersonic  Flows 

Virtually  throughout  the  shock  layer  the  flow  is  out  of  equilibrium.  However,  as  the 
velocity  approaches  zero  at  the  stagnation  point  the  associated  convective  time  scale  increases 
without  limit  Since  die  chemical  time  scale  remains  finite,  local  equilibrium  is  achieved  at  the 
stagnation  point  in  the  steady  state.  At  all  points  in  the  field  the  Law  of  Mass  Action 
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determines  a  local  equilibrium  concentration ,  ce,  which  "drives"  the  actual  concentration  c . 
Thus,  the  local  concentration  must  reach  the  local  equilibrium  concentration  at  the  stagnation 
point. 


In  terms  of  a  simple  Landau-Teller  type  model  for  the  chemical  source  term  W , 

W-Ce~.c 

tchtm  (6.2) 


FIGURE  6.7 

Schematic  of  generalized  concentration  and  source  term  behavior  along  the  stagnation  streamline  in  inviscid  now. 

In  Equation  (6.2)  the  degree  of  nonequilibrium  governs  the  magnitude  of  the  source  term. 
Figure  6.7  contains  a  schematic  of  the  stagnation  streamline  process  in  terms  of  the  c  and  W 
behavior. 

The  stagnation  point  is  at  the  right  of  both  plots  (x/8=  0.0).  Since  the  shock  itself  is 
frozen,  initially  a  large  difference  can  exist  between  c  and  ce.  In  region  1,  the  flow  is  far  out 
of  equilibrium  and  Tchem  is  very  small,  making#  large.  Thus,  the  concentration  adjusts 
rapidly  as  dissociation  soaks  up  internal  energy.  As  this  process  lowers  the  temperature,  the 
exponent  in  the  forward  reaction  term  (Eq.  2.33)  decreases  rapidly,  changing  W  in  region  2 
by  several  orders  of  magnitude.  As  tchem  increases,  #  decreases  at  approximately  the  same 
rate,  (Eq.  6.2),  and  the  progress  of  c  toward  its  goal  of  ce  is  slowed. 

The  rising  density  in  the  flow  approaching  the  stagnation  point  increases  the  number  of 
particle  collisions,  thus  elevating  the  importance  of  the  backward  rate  term  in  the  source  term 
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expression.  The  competition  between  backward  and  forward  terms  as  they  equalize  drives  W 
to  zero,  and  c  approaches  ce . 


FIGURE  6.8 

Stagnation  streamline  profiles  of  atom  mass  fraction  and  local  equilibrium  mass  fraction  for  the  blunted  wedge  detailed 
in  Table  6.1  with  chemical  rate  retarded  two  orders  of  magnitude. 


Figure  6.8  contains  a  computed  result  for  a  blunted  wedge.  The  behavior  along  the 
stagnation  streamline  and  its  continuation  along  the  body’s  surface  is  shown.  Conditions  for 
this  case  match  the  comparison  with  computation  detailed  in  Table  6.1,  but  here  the  chemical 
rate  was  retarded  by  two  orders  of  magnitude  to  emphasize  nonequilibrium  behavior. 

The  figure  illustrates  aspects  discussed  after  Fig.  6.7,  and  in  particular,  the  final  approach 
to  equilibrium  just  before  the  flow  stagnates  at  x/Rn  *  0.  The  species  equations  do  not  solve 
the  Law  of  Mass  Action  directly,  but  rather,  imply  a  steady  state  solution  when  the  source 
terms  vanish  at  the  stagnation  point .  Since  ce  at  the  stagnation  point  comes  from  a  direct 
solution  of  the  Law  of  Mass  Action,  it  is  an  independent  check  of  the  overall  validity  of  the 
solution  to  the  governing  equations.  A  slight  discrepancy  in  stagnation  point  concentration 
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results  from  the  presence  of  numerical  dissipation  and  disappears  in  the  limit  of  either  infinite 
resolution  or  zero  smoothing. 


Symmetry  plane  profile  of  the  logarithm  of  the  source  term  normalized  by  that  just  behind 
the  frozen  shock  for  the  blunted  wedge  detailed  in  Table  6.1  with  chemical  rate  retarded  two 
orders  of  magnitude. 

As  the  flow  expands  around  the  nose,  the  density  drops  rapidly,  reducing  the  collision 
rate  and  effectively  freezing  out  the  chemistry.  In  this  region,  W  remains  small  in  part 
because  contributions  from  both  the  forward  and  backward  terms  remain  relatively  small. 

Consistent  with  the  source  term  remaining  small  compared  with  its  post  shock  value,  the 
chemical  time  scale  remains  quite  large.  The  concentration  tends  toward  its  local  equilibrium 
value  very  slowly  and  would  achieve  it  only  far  downstream. 

Effects  sfJBisgwiatign  Energy 

As  the  symmetry  plane  flow  dissociates  downstream  of  the  normal  shock,  the  growth  of 
the  dissociation  energy,  a&o,  portion  of  the  internal  energy  implies  a  relatively  decreasing  T. 
The  appearance  of  &d  within  the  exponential  factor  of  the  source's  forward  reaction  suggests 
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that  for  some  net  reaction  rate  and  associated  chemical  length  scale,  the  So  role  is  crucial  to 
the  extent  of  chemical  length  scale  changes  throughout  the  shock  layer.  For  large  &o  the 
a&o  contribution  to  the  internal  energy  expression  may  dictate  a  large  temperature  change  for  a 
small  change  in  concentration.  The  source  term  will  then  "feel"  this  temperature  change  after 
being  selectively  amplified  through  the  exponential  in  the  forward  rate  term.  This  mechanism 
produces  radical  changes  in  chemical  time  scale,  Eq.  (6.2).  Moreover,  since  only  a  small 
change  in  c  produced  these  changes,  c  may  remain  far  out  of  equilibrium. 

Examine,  for  example,  Figure  6.9  which  shows  the  normalized  source  term  along  the 
symmetry  plane  on  a  logarithmic  scale  for  the  same  case  as  Figure  6.8.  Despite  the  fact  that 
over  the  first  five  points  downstream  of  the  shock  the  concentration  change  is  only  about  20% 
(  Fig  6.8),  the  source  term  changes  by  an  order  of  magnitude.  That  is,  rChem  changes  by 
approximately  an  order  of  magnitude.  This  implies  a  much  longer  relaxation  length,  and 
although  the  flow  remains  far  from  equilibrium,  the  progress  of  c  toward  ce  slows  radically. 

Of  course,  there  always  exists  some  overall  chemical  rate  or  Damkohler  number  large 
enough  to  eliminate  such  effects,  but  calculations  along  the  sustained  flight  corridor  suggest 
that  this  phenomenon  remains  important  -  especially  in  the  case  of  dissociating  nitrogen  at  low 
levels  of  dissociation.  Dissociating  N2  absorbs  nearly  twice  the  energy  of  dissociating  O2,  re¬ 
sulting  in  much  stronger  coupling  between  chemical  and  gas  dynamical  variables. 

As  c  slows  in  its  progress  toward  equilibrium  the  temperature  changes  more  slowly, 
decreasing  the  rate  at  which  W  changes  and  preventing  the  rapid  adjustment  of  c  toward  ce. 
Finally,  as  the  convective  time  scale  rises  in  the  stagnation  region,  the  finite  chemical  time  scale 
becomes  small  by  comparison,  producing  a  small  "boundary  layer"  region  of  rapid  adjustment 
as  the  concentration  equilibrates  with  ct.  This  effect  often  appears  as  the  "tail"  shown  in  region 
3  of  the  schematic  in  Figure  6.7,  and  again  displayed  in  the  computation  of  6.8.  Since  it  is  so 


109 


intimately  related  with  length  scale  behavior,  the  effect  is  seen  most  clearly  in  flows  with  large 
So  parameters. 

Effects  of  Reaction  Rate  on  Length  Scale  Behavior 

The  dissociation  rate  governs  the  rate  of  change  of  the  chemical  relaxation  length 
throughout  the  gas,  and  the  non-dimensional  reaction  rate  parameter  sets  the  magnitude  of  this 
scale.  The  previous  discussion  can  be  generalized  to  include  a  wide  class  of  nonequilibrium 
problems  ranging  from  nearly  frozen  to  nearly  equilibrium. 

The  sketch  in  Figure  6.10  shows  example  distributions  of  concentration  and 
corresponding  local  equilibrium  concentration  for  blunt  body  flows  with  Damkohler  numbers 
from  H*  0  (frozen)  to  -*  <»  (equilibrium  ).  If  the  flow  is  very  nearly  frozen,  ©,  the 
shock  layer  experiences  little  to  no  dissociation  upstream  of  the  stagnation  point,  where  it  must 
eventually  adjust  to  equilibrium.  Curves  ©  and  ©  exemplify  typical  profiles  found  in  flows 
for  f'of  approximately  0.1  and  10  respectively.  Numerical  calculations  at  these  conditions  are 
shown  later.  Finally,  ©  traces  the  concentration  behavior  of  flows  very  nearly  in  equilibrium. 
Here  the  infinitesimal  chemical  length  scale  permits  very  rapid  convergence  of  c  and  ce  and 
these  remain  equilibrated  until  the  flow  reaches  the  stagnation  point.  If  sufficiently  near  to 
identical  equilibrium,  these  distributions  will  remain  indistinguishable  throughout  the 
expansion.  Since  the  stagnation  region  is  normalized  by  the  stand-off  distance  6,  this  sketch 
does  not  show  how  the  shock  layer  thickness  increases  for  more  nearly  frozen  cases.  Over  the 
body  surface,  the  exact  local  equilibrium  concentration  will  depend  on  the  precise  values  of 
other  flow  variables.  However,  these  profiles  will  all  display  the  same  qualitative  behavior, 
and  the  sketch  presents  only  one  curve  downstream  of  the  stagnation  point. 
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equilibrium. 

Small  Departures  from  Equilibrium 

Figure  6.11  shows  the  behavior  of  the  wedge  example  presented  earlier  (Table  6.1  and 
Figure  6.3)  in  comparison  to  the  experiment  in  Reference  25.  Despite  the  near  equilibrium 
nature  of  the  flow  along  the  symmetry  plane,  the  length  sclIc  changes  enough  by  the  stagnation 
point  to  freeze  out  over  the  expansion  and  wedge  surface.  Notice  that  the  upward  climb  of  ce 
toward  c  results  from  an  increase  in  temperature  as  i7  atom.'  slowly  recombine  and  result  in  a 
slow  drop  in  c.  Capturing  the  initial  transient  in  a  simulation  with  this  severe  a  discrepancy 
between  chemical  and  convective  length  scales  requires  very  high  resolution  near  the  shock. 
The  example  illustrates  the  need  for  extreme  caution  before  awarding  identical  equilibrium 
procedures  to  flows  with  high  dissociation  energies.  Here  we  see  that  even  in  extremely  high 
Damktthler  number  flows,  the  chemical  length  scale  may  change  enough  to  freeze  out 
appreciably. 


Ill 


Stagnation  streamline  mass  fraction  and  local  equilibrium  for  circularly  blunted  IS*  wedge. 
(See  also  Table  6.1  and  Figure  63). 


Removal  of  Stagnation  Singularity  -  Viscous  Cases 

Curve  ®  in  6.10  demonstrates  near-singular  behavior  as  the  frozen  concentration 
through  the  shock  layer  rapidly  advances  toward  equilibrium  at  the  stagnation  point  Inviscid 
calculations  demonstrate  this  behavior  since  the  surface  streamtube  remains  quite  hot.  Realistic 
viscous  computations  and  real  flows,  however,  maintain  a  relatively  cool  thermal  boundary 
layer  immediately  adjacent  to  the  surface  which  removes  this  singularity.  As  atoms  re-combine 
in  the  cool  layer  just  removed  from  the  surface,  the  actual  concentration  and  local  equilibrium 
tend  toward  zero  dissociation.  Of  course,  wall  catalytic  effects  modify  this  behavior  at  elevated 
wall  temperatures.  The  present  examples  maintain  a  wall  temperature  of  1,500  K,  cold  enough 
to  keep  catalytic  effects  miniscule. 

Figure  6.12  details  this  behavior  for  two  cases  which  correspond  to  ©  and  ©  in  Figure 
6. 10.  In  the  figure  on  the  left,  *¥  -  10,  while  at  the  right  =  0.1.  Both  cases  correspond  to 
Mach  12  flight  in  the  standard  atmosphere  at  60  km  altitude,  and  the  viscous  examples  use  Re 
=  6500,  Pr  =  0.72,  and  Sc  =  0.5. 
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T)£,  Normalized  Dirunce  From  Stic.  11/5.  Nomulized  Disunce  From  Sue.  Point 

FIGURE  6.12 

Mass  fraction  and  local  equilibrium  mass  fraction  profiles  along  the  symmetry  plane  for  viscous  and  inviscid  flows  at 
60  km  altitude  in  STD  atmosphere.  V  -  10  (left),  i1  =  0.1  (right),  Re  *6500,  Pr  *  0.72,  Sc  =  0.5. 

It  is  interesting  to  note  that  the  local  equilibrium  concentration  does  not  lead  the  actual  concen¬ 
tration  to  zero  in  the  stagnation  region  for  the  viscous  examples.  Behavior  in  this  region  stems 
from  heat  conduction  and  other  viscous  effects.  The  viscous  terms  on  the  right  side  of  the 
governing  equations  (Eq.  1.9)  act  as  additional  sources  and  the  simple  source  term  model  in 
Equation  (6.1)  neglects  these  complex  terms  altogether.  In  fact,  non-zero  behavior  stands  as 
evidence  that  these  terms  contribute  in  this  region.  Moreover,  the  fact  that  the  viscous  profiles 
trace  the  inviscid  curves  almost  precisely  through  the  inviscid  portion  of  the  field  suggests  that 
mass  diffusion,  heat  transfer,  and  shear  terms  remain  comparatively  small  through  this  portion 
of  the  field. 

For  purposes  of  comparison,  the  horizontal  scales  in  Fig.  6.12  were  normalized  by 
stand-off  distance.  This  normalization  masks  the  fact  that  the  thick  boundary  layer  in  the 
viscous  cases  actually  displaced  the  shock  by  ~7%. 
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6.3  Coupled  Reacting  Systems 


When  multiple  reactions  occur  in  a  gas  mixture,  the  problem  of  a  disparate  chemical  and 
convective  length  scale  becomes  one  of  many  disparate  length  scales.  The  degree  of  coupling 
between  reactions  in  such  a  system  depends  largely  on  two  factors.  First,  the  relative  amounts 
of  energy  absorbed  or  produced  by  a  reaction  contribute  directly  to  the  overall  internal  energy 
of  the  system  and  directly  affect  thermodynamic  properties.  For  example,  as  one  reaction 
proceeds,  it  may  absorb  so  much  internal  energy  that  some  other  reactions  freeze  out. 
Secondly,  the  relative  scales  of  the  processes  affect  reaction  coupling.  If  one  reaction  proceeds 
much  more  rapidly  than  another,  the  fast  reaction  may  reach  its  equilibrium  concentration  long 
before  the  other  reaction  has  begun  to  affect  the  internal  energy  of  the  gas. 


Formation  of  NO  wjtfiiiuhs  Shock  Layer 

In  reacting  air  systems,  such  coupling  effects  are  very  clearly  demonstrated  upon 
examination  of  the  formation  of  NO  throughout  the  shock  layer.  As  an  example,  re-consider 
the  axisymmetric  multi-reaction  test  case  presented  in  Section  6.1  (Table  6.2,  Figs.  6.4  and 
6.5).  Since  the  shock  layer  is  in  vibrational  nonequilibrium,  the  detailed  concentration  levels  of 
mixture  components  are  incorrect  by  the  standard  of  more  accurate  chemical  modeling. 
Nevertheless,  the  species  behavior  in  this  comparison  agrees  reasonably  well. 


Figure  6.13  provides  profiles  of  N,  O,  NO,  and  O2  along  the  stagnation  streamline.  At 
conditions  just  behind  the  shock  the  characteristic  relaxation  length  for  molecular  oxygen 
dissociation  is  about  0.01fl„,  while  that  for  nitrogen  is  approximately  0.5/?„.  This  50-fold 
disparity  is  a  measure  of  the  degree  of  coupling  between  the  two  dominant  reactions. 


114 


Mass  Fraction 


Symmetry  plane  profiles  of  species  mass  fraction  for  Oi,  0,  NO,  and  N  for  the  case  detailed  in  Table  6.2  (also 
see  Figs.  6.5  A  6.6). 

The  earlier  plot  of  equilibrium  constants  in  Figure  1.3  provided  an  indication  of  the 
relative  importance  of  the  forward  and  backward  terms  in  the  reaction  rate  expressions  for  the 
five  reactions  considered  here.  Of  course,  in  view  of  the  actual  chemical  nonequilibrium.  Fig. 
1.3  only  provides  an  approximate  guide  to  the  source  term's  behavior.  As  the  flow  in  Fig. 
6.13  crosses  the  Mach  13  normal  shock,  oxygen  rapidly  dissociates,  creating  an  abundance  of 
free  oxygen  atoms.  These  atoms,  in  tum,  activate  the  exchange  reactions  which  proceed  at 
rates  slightly  faster  than  that  of  O  production.  However,  the  competition  between  these 
reactions,  mentioned  in  the  closing  paragraphs  of  Chapter  1,  slows  the  net  production  of  NO  to 
about  Vioth  that  of  O2  dissociation.  At  these  post-shock  temperatures,  the  first  exchange  reac- 
tion  (O  +  N2  <=>  N  +  NO)  proceeds  faster  than  the  NO,  robbing  the  exchange  reaction  (NO  +  O 
<=>  N  +  O2).  It  is  this  mechanism  which  is  responsible  for  the  initial  production  of  NO  behind 
the  shock. 

As  the  temperature  decreases  to  -8,000  K  inside  the  shock  layer,  this  competition  results 
in  almost  exact  cancellation  of  the  NO  produced  and  destroyed  by  the  shuffle  reactions.  Thus 
the  only  mechanism  left  for  altering  NO  concentration  is  the  NO  dissociation  reaction.  At 
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these  temperatures,  the  forward  rate  term  dominates,  breaking  up  NO  molecules  at  a  rate 
slightly  faster  than  that  of  nitrogen  dissociation.  This  effect  results  in  the  decline  of  NO  con¬ 
centration  approaching  the  stagnation  point 

The  remainder  of  the  shock  layer  is  shown  in  Figure  6.14  for  mass  fraction  contours  of 


FIGURE  6.14 

Contours  of  NO  mass  fraction,  and  temperature  ratio  computed  for  the  case  detailed  in  Table  6.2 
(also  see  Figs.  6.5  &  6.6). 


The  temperature  field  holds  two  important  functions  in  this  discussion.  First  for  any  point  in 
the  field,  it  determines  values  of  the  equilibrium  constants  for  the  reactions,  giving  an  indica¬ 
tion  of  what  local  equilibrium  condition  the  species  concentrations  seek.  Secondly,  the  temper¬ 
ature  provides  an  indication  of  the  overall  magnitude  of  the  source  terms,  and  therefore  the 
relevant  chemical  length  scale  at  any  point  in  the  field. 

As  the  flow  expands  out  of  the  hot  stagnation  region,  the  temperature  decreases  rapidly, 
slowing  the  progress  of  all  reactions  toward  their  respective  equilibrium  condi tions*.Since 
nitrogen  dissociation  has  the  slowest  chemical  rate,  it  is  the  first  to  freeze  out.  Immediately 
afterward,  the  nitric  oxide  dissociation  reaction  follows  suit  With  the  dissociative  reaction 
turned  off,  the  exchange  reactions  exclusively  determine  the  concentration  of  NO  in  the  shock 
layer,  and  since  their  individual  rates  are  so  high,  these  reactions  continue  to  be  active,  even  in 
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the  rapidly  freezing  flow  over  the  body's  shoulder.  The  unbalance  of  these  reactions  continues 
to  drive  the  concentration  of  NO  upward  through  most  of  the  expansion  fan.  This  interesting 
point  also  was  documented  by  Park  (30),  and  is  a  direct  result  of  the  fact  that  the  first 
exchange  reaction  (O  +  N2  <=>  N  +  NO)  still  produces  more  NO  than  the  other  consumes. 
However,  as  the  temperature  continues  to  drop,  this  situation  reverses  itself,  resulting  not  only 
in  a  decrease  in  NO ,  but  also  O2  formation  (by  the  same  reaction).  As  oxygen  molecules 
continue  to  deplete  the  supply  of  free  oxygen  atoms,  NO  does  not  have  the  opportunity  to  re¬ 
form.  Of  course,  all  of  these  processes  take  place  over  much  expanded  lengths,  since  on 
reaching  the  shoulder  of  the  body  the  expanding,  cooling  flow  has  changed  the  chemical  length 
and  time  scales  by  several  orders  of  magnitude. 
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7 .  Effectiveness  of  Adaptive  Grid  Embedding  in 
Hypersonic  Flows 


The  test  cases  and  investigation  demonstrate  the  feasibility  of  applying  the  adaptive 
methodology  to  hypersonic  reacting  flows.  In  evaluating  the  effectiveness  of  this  procedure, 
several  aspects  require  consideration.  A  comparison  of  such  mechanical  points  as  computation 
time  and  memory  savings  afforded  by  the  present  technique  are  of  interest,  as  are  the  resolution 
requirements  of  some  of  the  physical  phenomena  captured  in  embedded  regions  of  the  flow 
field  and  possible  improvements  in  the  implementation  of  the  technique  itself. 


7.1  Enhancement  of  Computational  Efficiency 

Since  the  purpose  of  the  final  grid  is  to  resolve  all  length  scales  in  the  domain,  the 
computational  savings  associated  with  any  flow  field  depend  heavily  on  the  number  of  length 
scales  involved.  Figure  7.1  contains  a  highly  adapted  grid  resulting  from  a  Mach  16,  2-D 
simulation  at  conditions  45  km  aloft  in  the  standard  atmosphere.  The  conditions  were  modified 
to  enhance  oxygen  dissociation  by  artificially  freezing  the  nitrogen  and  nitric  oxide  dissociation 
reactions  {Osi  =  &NO  a  0)-  The  body  size  is  such  that  the  characteristic  length  for  oxygen 
dissociation  just  downstream  of  the  normal  shock  is  only  0.001  Rn  and  virtually  all  of  the 
molecules  dissociate  within  1%  of  the  nose  radius. 
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FIGURE  7.1 

Adapted  grid  with  four  levels  of  embedding  and  1550  nodes  for  conditions  described  in  text. 


Although  somewhat  contrived,  this  case  provides  a  good  demonstration  of  the  capabilities 
of  the  embedded  mesh  procedure.  Moreover,  examining  the  time  savings  after  each  level  of 
embedding  provides  a  realistic  estimate  of  the  savings  for  a  wide  range  of  chemical  rates. 

Table  7. 1  lists  the  computation  times  required  for  the  original  coarse  grid,  final  adapted 
grid,  and  four  globally  embedded  grids.  Notice  that  at  the  two  highest  levels  of  global  mesh 
refinement,  estimated  simulation  times  became  sufficiently  long  that  actual  tests  were 
precluded. 

The  final  adaptive  grid  resolved  the  chemical  relaxation  with  four  levels  of  embedding, 
and  the  computation  cost  proved  to  be  8.6  times  that  of  the  coarse  grid  baseline  case. 
Resolving  the  rapid  decay  on  a  globally  fine  mesh  would  require  the  same  four  embeddings, 
but  at  a  computational  cost  exceeding  250  times  the  baseline  computation. 
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TABLE  7.1 


Computing  time  comparisons  for  adaptive  and  globally  refined  meshes  for  conditions  of  Figure  7.1. 


Mesh 

Dimensions 

#of 

nodes 

Normalized 
Computer  Time 

Original  Coarse  Grid 

10x20 

200 

1.0 

Adapted  (4-level) 

Unstructured 

1550 

8.6 

Entire  mesh  (1  level  embedded) 

19x39 

741 

4.7 

Entire  mesh  (2  levels  embedded) 

2849 

17.9 

Entire  mesh  (3  levels  embedded) 

11322 

70.8* 

Entire  mesh  (4  levels  embedded) 

43617 

274.7* 

♦Estimated  due  to  computing  resource  limitations 


Such  comparisons,  however,  are  somewhat  optimistic.  Obviously,  a  globally  fine  grid 
would  provide  exceptional  detail  throughout  the  field,  while  the  adaptive  solution  only  provides 
resolution  comparable  to  the  scale  of  the  flow  features  in  particular  regions.  Moreover,  given 
an  initial  solution  on  a  coarse  grid  to  identify  flow  features,  one  could  conceivably  apply  more 
traditional  grid  control  on  a  structured  mesh  providing  adequate  resolution  of  the  relaxation 
layer  with  far  fewer  than  the  44,000  nodes  required  by  the  globally  refined  mesh.  The  net  time 
for  computing  a  solution  with  such  a  procedure,  however,  is  not  clear.  The  final  grid  would 
not  place  nodes  as  accurately  as  the  adaptive  grid,  and  would  certainly  contain  more  nodes. 
The  time  required  to  set  up  such  a  grid  is  difficult  to  estimate  quantitatively.  Finally,  since  such 
a  comparison  is  so  intimately  related  to  the  chemical  relaxation  length,  it  would  be  all  but 
meaningless  as  an  indication  of  typical  adaptive  savings. 

As  a  better  measure  of  the  savings  associated  with  adaptation,  consider  the  examples  put 
forth  earlier.  Reference  25  computed  the  blunt  wedge  case  on  3000  nodes  to  provide  the 
resolution  shown  in  Figure  6.3.  However,  that  figure  makes  clear  that  the  current  adaptive 
technique  provides  superior  resolution  with  just  over  half  the  computational  nodes.  The  coarse 
base  grid  for  this  case  contained  300  nodes  and  was  chosen  to  provide  adequate  resolution  of 
frozen  flow  features.  The  adaptive  solution  required  5.7  times  the  computation  effort  of  the 
coarse  grid  solution  to  converge  (5  orders  of  magnitude  in  RMS  momentum  residual). 
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The  axisymmetric  sphere  results,  discussed  extensively  in  Section  6.3  (Figs.  6.12, 
6. 13),  present  further  evidence  of  the  cost  associated  with  adaptive  hypersonic  solutions.  The 
final  grid  for  this  case  contained  1300  nodes,  and  required  three  levels  of  embedding  and  S.3 
times  the  computing  time  of  the  original  200  node  mesh. 

These  results  typify  the  experience  gained  over  the  course  of  this  investigation.  Shock 
fit,  reacting  hypersonic  flow  fields  over  simple  geometries  usually  require  5-6  times  the 
computational  effort  of  a  coarse  base  grid  solution.  In  this  context,  a  "coarse  grid"  is  one 
which  clearly  shows  frozen  flow  features,  but  whose  resolution  is  only  about  half  of  that 
required  for  engineering  computations.  The  coarse  base  grids  referred  to  in  the  preceding  two 
paragraphs  exemplify  these  qualities.  As  the  example  in  Figure  7.1  suggests,  flows  with  larger 
disparity  between  chemical  and  body  length  scales  typically  require  more  computational  time, 
while  those  more  nearly  fro 7-  quire  less.  Remember  that  frozen  calculations  to  engineering 
precision  take  about  4  tun-  more  effort  to  compute  than  the  coarse  base  grid  referred  to  by  this 
comparison  (see  the  "entire  mesh  [1  level  embedded]"  entry  in  Table  7.1).  Thus,  the  time 
required  for  a:,  adaptive  blunt  body  reacting  flow  field  is  approximately  1.5  times  that  needed 
for  a  frozen  solution  (using  globally  refined  meshes)  of  comparable  accuracy. 

This  is  an  extremely  important  point,  with  respect  to  the  effectiveness  of  the  adaptive 
method.  The  example  shown  in  Figure  7.2  makes  this  point  even  clearer.  At  the  left,  the  fig¬ 
ure  shows  temperature  ratio  contours  for  inviscid  Mach  12  frozen  flow  over  an  axisymmetric 
body.  The  grid  required  800  nodes  to  adequately  resolve  the  frozen,  inviscid  features,  and 
since  frozen,  contains  only  one  length  scale. 

The  right  of  Figure  7.2  shows  an  adapted  case  computed  in  uncoupled  (NO  absent) 
reacting  air  at  the  same  Mach  number.  At  30  km  altitude  in  the  standard  atmosphere,  a  nose 
radius  of  0.0318  m  produces  a  characteristic  chemical  relaxation  length  of  approximately 
0.01  Rn  behind  the  normal  shock.  This  second  length  scale  greatly  increases  the  resolution 
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requirements  of  the  numerical  procedure.  Nevertheless,  two  levels  of  adaptive  embedding 
produced  results  with  resolution  comparable  to  that  of  the  frozen  case.  The  final  adapted  grid 
contains  1341  nodes  and  resolves  both  the  chemical  and  body  length  scales.  Here  the  adaptive, 
nonequilibrium  solution  required  a  1.6  times  greater  computation  time  than  a  frozen  simulation 
over  the  same  geometry  -  despite  the  multiple  length  scales  in  the  real  gas  case. 


FIGURE  7.2 

Contoun  of  temperature  ratio  for  inviicid  Mach  12  flow  for  frozen  flow  (left),  and  meeting  flow  (right). 


7.2  Adaptive  Resolution  of  Physical  Phenomena 

Most  algorithms  capture  shocks  in  the  interior  of  a  computational  domain.  This  process 
usually  results  in  a  smearing  out  of  the  physical  disturbances  over  several  nodes.  Consider, 
for  example,  the  shock  captured  solutions  of  Reference  5  presented  in  the  first  section  of 
Chapter  6  (Figs.  6.5,  6.6).  There,  shock  capturing  has  rounded  off  the  density  and 
temperature  jumps  across  the  bow  shock.  Obviously,  an  increase  in  grid  resolution  near  the 
discontinuity  would  result  in  sharper  shock  jumps,  but  without  a  priori  knowledge  of  shock 
positions,  this  would  be  expected  to  be  an  expensive  process. 
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Shock  Triggered  Noneouilibrium 


Imagine  a  flow  field  with  a  captured  shock  in  a  hypersonic  flow.  Since  the  calculation 
may  need  to  be  performed  at  several  angles  of  attack,  and  since  chemical  effects  will  alter  shock 

i 

shape,  exact  shock  positions  may  not  be  found  prior  to  computations.  However,  at  hypersonic 
stagnation  enthalpies,  any  shocks  may  trigger  chemical  reactions.  Since  these  reactions 
determine  the  downstream  gas  composition,  the  relaxation  layer  requires  adequate  resolution. 

Now,  re-examine  the  shock  captured  density  and  temperature  profiles  in  Figures  6.5  and 
6.6a.  The  relaxation  zone  is  an  order  of  magnitude  wider  than  the  shock  smearing,  but  for  ar¬ 
bitrary  shock  conditions,  this  will  not  always  be  the  case.  Consider  a  scramjet  engine  inlet,  in 
which  each  captured,  reflected  shock  will  increase  the  static  temperature,  decreasing  the 
relaxation  length  until  the  relaxation  zone  may  reside  entirely  within  the  captured  shock.  Such 
a  situation  would  obviously  be  detrimental  to  accurate  predictions  of  downstream 
concentrations  and  flow  properties.  References  33  and  34  contain  an  examination  of  adaptive 
solutions  through  such  geometries  and  demonstrate  the  adaptive  scheme's  ability  to  separate  the 
chemical  relaxation  from  the  translational  shock. 

A  second  example  of  such  shock-induced  nonequilibrium  effects  occurs  in  trans- 
atmospheric  or  re-entry  flight.  In  this  regime,  the  translational  temperature  behind  the  bow 
shock  will  be  high  enough  to  incite  ionization  before  dissociation  and  vibration  are  able  to 
decrease  this  extreme  temperature.  Accurate  prediction  of  the  electron  density  in  the  gas  cap 
therefore  depends  wholly  upon  the  temperature  spike  resolved  after  the  translational  shock  but 
before  nonequilibrium  thermo-chemical  processes  erode  this  temperature.  Again,  the  multiple 
length  scales  of  these  processes  ideally  suit  such  problems  to  adaptive  computations. 
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7.3  Recommendations  for  Improving  the  Adaptive  Technique 

The  examples  and  test  cases  in  Chapter  6  demonstrate  the  adaptive  scheme's  ability  to 
locate  and  resolve  flow  features.  In  the  current  implementation,  the  codes'  suggested 
adaptation  patterns  almost  always  enclose  the  detected  feature,  placing  cell  interfaces  in 
relatively  benign  portions  of  the  flow  outside  the  critical  area.  This  section  suggests 
enhancements  to  the  basic  adaptation  methodology  based  on  experience  with  the  existing 
codes. 

Directional  Embedding 

Figure  7.3  displays  two  examples  of  directional  embedding.  Even  in  complex 
hypersonic  flow  fields,  flow  features  often  align  with  the  body  and  coarse  grid  orientation. 
Directional  embedding  takes  advantage  of  this  alignment  to  increase  solution  resolution  while 
still  avoiding  the  creation  of  unnecessary  cells  and  nodes.  For  example,  upon  examination  of 
the  adapted  grids  presented  in  earlier  chapters,  it  becomes  clear  that  many  of  the  features 
discussed  are  largely  one  dimensional.  Both  the  chemical  relaxation  zone  and  boundary  layer 
remain  primarily  grid  aligned  and  would  benefit  from  such  directional  adaptation.  Kallinderis 
(17)  investigated  this  subject  in  some  detail  and  claims  a  reduction  in  computational  effort  of 
approximately  two  (over  an  adaptive  solution  without  directional  embedding)  as  a  result  of  the 
fewer  cells  and  nodes  created  with  this  method  of  cell  division. 

Incorporation  of  directional  embedding  into  the  current  algorithm  requires  relatively  few 
changes.  The  feature  detection  algorithm  would  then  store  differences  in  both  directions  for 
each  cell.  Then  the  algorithm  for  cell  creation  and  pointer  updating  must  reflect  the 
unidirectional  nature  of  the  embedding. 
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▼  w - j—  — r - r 

FIGURE  7J 

Cell  subdivisions  for  directional  embedding. 


Adaption  Parameters 

Figure  7.4  shows  a  typical  adapted  grid  for  a  2-D  reacting  viscous  calculation  (Re  = 
6500).  This  grid  shows  embedded  regions  primarily  near  the  shock  and  body,  capturing  both 
the  relaxation  layer  and  boundary  layer  respectively. 


FIGURE  7.4 

Adapted  grid  for  a  2-D  reacting  viscous  calculation  and  total  velocity  vectors  for  same  calculation. 

The  current  use  of  density  and  concentration  differences  as  adaptive  parameters  stems 
from  expectations  of  inviscid  flow  field  features.  In  viscous  hypersonic  problems  however, 
the  thermal  gradients  near  the  wall  give  rise  to  density  gradients  via  the  equation  of  state. 
Thus,  density  differences  retain  their  utility  in  detecting  such  boundary  layer  behavior. 
However  it  is  the  thermal  boundary  layer  that  these  differences  indicate  (see  velocity  field  in 
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Fig.  7.4),  and  clearly,  density  differences  are  not  appropriate  for  general  shear  layers. 
Reference  17  has  made  use  of  both  velocity  gradients  or  shear  stresses  as  logical  adaptive 
parameters  to  identify  such  a  feature. 
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Summary  and  Conclusions 


The  application  of  adaptive  grid  embedding  to  hypersonic  flows  was  motivated  by  the 
scales  associated  with  nonequilibrium,  high  temperature  gases  that  must  be  modeled  for  such 
fields.  The  adaptive  technique  permitted  numerical  simulation  of  nonequilibrium  chemically 
reacting  flows  with  only  1.5  to  2  times  the  computational  effort  of  comparably  resolved  frozen 
flows.  Example  calculations  were  completed  for  two  dimensional,  inviscid,  flows  and 
extensions  to  both  axisymmetric  and  viscous  flow.  In  all  cases  the  algorithms  provided  a 
means  of  computing  and  analyzing  complex,  reacting,  hypervelocity  flows  over  simple 
geometries  such  as  blunted  wedges  and  cones. 

The  primary  contributions  of  this  work  relate  to  both  the  development  of  the  numerical 
algorithm  and  some  details  of  physical  phenomena  that  adaptive  gridding  makes  clearer. 

Numerics  and  Algorithm  Development 

This  work  extended  the  explicit  Lax-Wendroff  technique  developed  by  Ni  (29)  to 
hypersonic  calculations  on  an  unstructured  grid.  The  axisymmetric  computational  results  are 
believed  to  demonstrate  the  first  extension  of  Ni's  two-dimensional  algorithm  to  circular 
cylindrical  coordinates. 

In  addition  to  enhancements  of  the  integration  scheme,  the  nonequilibrium  shock-fit 
domain  necessitated  several  major  modifications  of  the  original  adaptive  method  (6).  For 
example,  shock  fitting  required  re-mapping  of  the  unstructured  computational  grid  at  each  time 
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step.  The  simple,  effective  technique  presented  in  Section  5.3  accomplishes  this  task  on 
unstructured  domains  with  any  number  of  embedded  levels  while  preserving  conservation 
across  mesh  interfaces.  A  lesser  contribution  to  the  adaptive  methodology  was  the  introduction 
of  a  feature  detection  algorithm  capable  of  triggering  mesh  refinement  dependent  upon  an 
arbitrary  number  of  nonequilibrium  and  flow  variables. 

Gas  Dynamics 

A  mixture  gas  model  was  developed  based  on  the  assumption  of  a  half-excited  vibrational 
state  for  diatomic  molecules.  This  model  included  N2, 02,  NO,  N,  and  O,  and  its  development 
required  the  identification  of  a  characteristic  density  for  dissociation  of  NO  molecules.  Section 
2.2  detailed  the  formulation  of  this  parameter.  Figure  2.2  supports  the  hypothesis  of  Pdno  ~ 
const,  required  by  the  dissociative  reaction  model. 

Both  the  two-dimensional  and  axisymmetric  flow  solvers  were  validated  through  perfect 
gas,  dissociating  gas  and  multiple  reaction  comparisons  and  an  investigation  of  hypersonic 
flow  over  blunted  cones  and  wedges  for  the  spectrum  of  reactions  varying  from  near  frozen  to 
near  equilibrium.  Examination  of  chemical  and  convective  time  and  length  scales  in  the 
domain,  provided  insight  into  the  character  of  inviscid  and  viscous  reacting  flows  under  widely 
varying  conditions  of  nonequilibrium. 

Of  special  note  is  the  interplay  between  dissociative  and  exchange  reactions  in 
nonequilibrium,  high  temperature  air  and  the  role  of  the  chemical  length  scale  in  reaction 
coupling.  For  example,  increases  in  the  concentration  of  NO  off  the  symmetry  plane  are 
linked  directly  to  the  effects  of  length  scale  and  reaction  coupling.  Finally,  some  viscous 
investigations  demonstrated  that  these  same  concepts  retained  their  importance  in  viscous 
hypersonic  flows.  However,  in  regions  dominated  by  diffusion  processes,  the  local 
equilibrium  concentration  must  share  its  "driver"  role. 
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Conclusions 


Unstructured,  adaptive,  embedded  gridding  has  been  applied  to  hypersonic,  blunt  body, 
nonequilibrium,  CFD  problems.  The  high  temperature  gas  mixture  was  described  by 
Lighthill's  dissociating  gas  model  which  was  extended  to  include  five  species  and  multiply 
coupled  reaction  paths  in  both  viscous  and  inviscid  flows.  The  Lax-WendrofF  numerical  algo¬ 
rithm  was  extended  to  include  shock  fitting  and  adaptation  on  general,  unstructured,  moving 
(adjusting)  grids. 

Comparisons  with  experiment  and  computation  provided  validation  for  both  the  explicit 
real-gas  algorithm  and  unstructured  shock  fitting  procedure  for  perfect  gas,  dissociating  gas, 
and  multiple  reaction  cases.  The  dissociating  gas  model  demonstrated  good  agreement  with 
shock  tunnel  experiment  and  computational  results  of  Reference  19,  and  reasonable  prediction 
of  gas  dynamic  variables  in  flows  with  coupled  reactions. 

A  detailed  study  of  basic  nonequilibrium  flow  phenomena  has  been  completed  for 
freestream  Mach  numbers  from  5  to  15  over  blunted  cones  and  wedges.  These  flows 
demonstrated  degrees  of  nonequilibrium  ranging  from  nearly  frozen  to  near  equilibrium. 
Adaptation  proved  useful  in  examining  behavior  along  the  stagnation  streamline,  especially  in 
cases  displaying  a  small  departure  from  equilibrium.  Here,  adaptation  was  shown  to  be 
particularly  useful  in  capturing  the  steep  chemical  gradients  which  appear  within  the  shock 
layer. 

Grid  adaptation  appears  to  be  a  cost-effective  way  of  computing  high  resolution  solutions 
to  hypersonic,  finite  rate,  real-gas  problems.  The  computational  effort  required  for  an  adap¬ 
tively  refined  nonequilibrium  solution  was  shown  to  be  only  1.5  to  2  times  that  required  for  an 
equivalently  resolved  frozen  flow  solution  for  the  same  configuration.  This  compares  favor- 
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ably  with  the  40-50  times  that  may  be  required  for  a  solution  on  a  globally  refined  mesh 
capable  of  resolving  the  chemical  relaxation. 
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APPENDIX  A 


Integration  Formulae 

A.1  Non- Orthogonal  Two-Dimensional  Coordinates 


Cell  Change  (refer  to  Fig  3.3) 


At 


AUc  =  AfWc-f=- 
Ac 


’+  -  x,) 

+  -  y.)-(GjtG|)(xj  -  x,) 


+ (ELri)0’'  -  »>  - -  *>J 


Ac  =  -  Xifoj  -  yi)-(Xj  -  xfok  -  yfl 


(A.l) 

(A.2) 


Distribution  Formulae  (refer  to  Fig.  3.4) 

8^4  =  ilAUc  “  Afc  ~  Age  +  yAWc 
5u4  =  ^|aUc  +  Afc  -  Age  +  yAWc] 
5Uc)t  =  j(AU  c  +  Afc  +  Age  +  yAWc] 
SUc),  =  J[aUc  -  Afc  +  Age  +  yAWc] 


where 
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Afc  —  ”t — (AFcAy^  —  AGcAxO 
Ac 

Age  =  — — (AGcAxm  —  AFcAym) 
Ac 


(A.4) 


Ax'  =  +  x*  -  x,  +  xj) 

&yl  =  fyyi  +  yt-  y»  +  yj) 

Atm  =  1<X*  +  Xy  -  X/  +  X.) 

Aym  =  £{y*  +  yy  -y/  +  y;) 


(A.5) 


Time  Step 


At  =  [CFL)  min||- 


_ Ac  _ k _ 

uAyt  —  vAx1  +  aAl\'  I uAym  —  vAxm  +  aAm\ 


(A.  6) 


In  A.6,  0  SCFL<,  1. 


Additionally, 


A /  =  V(Ax')2  +  (Ay')2,  Am  =  V  (Ax'")2  +  (Ay"*)2 


(A.7) 
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A.2  Non- Orthogonal  Axisjmmetric  Coordinates 


FIGURE  A.1 

A  general  cell  in  circular  cylindrical  coordinates. 


Referring  to  the  above  sketch  of  a  general  cell  in  an  axisymmetric  domain,  the  in  viscid  gov¬ 
erning  equations  in  circular  cylindrical  coordinates  become: 

aU+3F+i3[rG) 

d'  dx  '  dr  "  (A.8) 

where 
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Cell  Change 


where, 

r  =  ^r,  +  ry  +  r*  +  r/) 

Distribution  Formulae 

8Uc)i  =  *[aUc  -  AfC)  -  Agc,+  yAWc 
SUc|  =  J[aUc  +  At Cj  -  Age,  +  yAWc 
Slid,  =  J[aUc  +  Afct  +  Age,  +  yAWc 
8UC),  =  ijAUc  -  Afc,  +  Age,  +  yAWc 

with 

Afc„  -  ^(AFcAy/  -  jrAGc&x1)  .and  n  =  (ij,k,l) 
Ag cn  =  ^-^AGcAxm  -  AFcAym) 


(A.9) 


(A.  10) 


(A.  11) 


(A.  12) 


Along  the  axis  of  symmetry,  the  governing  equations  become 


The  cell  change  for  cells  along  this  axis  has  no  contribution  from  the  face  on  the  axis,  and  this 
equation  leads  to  Agc„  =  0  when  n  lies  along  axis  of  symmetry.  These  two  changes  com¬ 
pletely  account  for  the  singularity  in  the  coordinate  system  at  r  =  0. 

A.3  Viscous  Formulation  in  Two-Dimensions 


Primary  Cell  Secondary  Cell 


FIGURE  A.2 

Primary  and  Secondary  cells  for  viscous  integration. 

After  first  distributing  the  inviscid  changes  to  each  node  (via  Eq.A.3),  the  algorithm 
computes  nodal  changes  from  the  viscous  terms  R  and  S  in  the  governing  equations  (Eq.1.9). 
The  formulation  of  these  changes  is  naturally  more  complex,  since  these  dissipadve  terms  rely 
upon  second  derivatives.  As  with  the  inviscid  changes,  contributions  from  all  cells 
surrounding  a  node  complete  the  difference  stencil  for  these  terms.  The  Lax-Wendroff  scheme 
relies  upon  the  second  order  changes  of  the  inviscid  integration  for  stability.  Fortunately  this 
requirement  does  not  extend  to  the  viscous  terms  and  no  Jacobians  of  R  or  S  need  be 
computed  (18).  However,  this  omission  reduces  the  time  accuracy  of  the  scheme  to  first  order. 
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Following  the  notation  in  Figure  A.2  the  viscous  change  to  any  node  in  the  computational  mesh 
is 


+  [+  Ro&yB  +  RaW]-[+  SB&xB  +  S^Axfl 

+  [+  Rq  Ay£  +  ]  -  [+  S&  Ax£  +  ] 

+  [-  Rl^  -  RcW  3  -[-  SeBAxf  -  S£Ax?] 
+  [-*5Ari  -RM}-[-S2Ax1a  -  SgAxjj ] 


(A.  14) 


In  (1. 14),  the  superscripts  on  R  and  S  refer  to  cell  faces  and  the  area,  A,  is  the  area  of  the  pri¬ 
mary  cell  in  Figure  A.2. 

In  the  integration  scheme,  however,  we  sweep  through  the  cells  and  determine 
contributions  to  each  node.  Thus  Equation  A.  14  needs  to  be  re-expressed  in  terms  of 
contributions  from  any  cell  to  its  four  comer  nodes. 


Viscous  distribution  formulae 

(A  =  —[(+  RsAyl  -  R^Ay")  -  (+  S'  Ax'  -  S"  Ax'")] 

(At/vi,cU  *  7^(+  Rn  Ay*  +  Rw&ym)-(+  S'*  Ax'  +  SwAxm)] 
(A  C/vi,cW  =  Rn&yl  +  R*Aym)-(-  S"  Ax'  +  S' Ax'")] 

(A£/v«c)m  =  R'  Ay'  -  Re  Ay”)  -  (-  S'  Ax'  -  S' Ax'")] 

rl  • 


The  terms  in  the  R  and  S  vectors  in  this  equation  contain  first  derivative  terms  like  (uf)cd  which 


must  be  split  into  contributions  from  the  cells  involved. 


+  -^AUo'iyc  -  yo)  +  Uciybc  -  yc)  +  Uoiyo  -  y^B 
+  -^{Uo'iyo'  -  yo)  +  U0{y<ta  -  yo)  +  UD{yD-  -  y^)] 


(A.  16) 
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This  formulation  stresses  the  contributions  from  each  cell  C  and  D  in  the  composition  of  the 
first  derivatives. 

Reference  18  contains  a  thorough  explanation  of  the  viscous  integration  scheme  for  non- 
reacting  flow  in  general  two  dimensional  coordinates. 
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B.2  Noneqidlibrium  Mixture  Gas  Model 


For  the  nonequilibrium  mixture  gas  model  developed  by  Chapter  2  the  pressure  at  any 
point  in  the  field  is 

Where  Un  represents  the  /Ith  element  of  the  state  vector  and  jf  and  0  are  defined  by 

dtl  .  2$t\  +  ma) +  'Hwa +  ^4  + 


+  -  Victim,  +  Urj±  +  Us) 

-  UlCfracimz  +  Uljjfc  + 


with 


♦ 


C^acl  />7ti  +/2W2 
Cfi'ac2~f\mi  +hmi 


(B.7) 


4 

and  for  air/i  =  0.79  and /2  =  0.21. 
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Where 


^-  =  (1  -  CfracxmJ^  +  il  -  Cpaclm1)^  +  ^ 

d&  _  0pl 
dU5  ~  m* 

dS  _ 
dUj  ~  m* 

d&  -  l(^)i  .  _  &D} 

dlh~lSm3  »*3  I  m 


Jacobian  Matrices 
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0 

0 

0 

0 

0  “ 

dp  u\ 

2^+i2_ 

dp 

dp 

dp 

dp 

dp 

dV\  t/f 

u  1  ac/2 

ac/3 

ac/4 

dUs 

dU6 

dih 

Mi 

Ui 

0 
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0 

0 

U\ 

Ui 

t/i 

u*m 

Uim+h° 

{/2 1^2. 
Wj 

< 

C/2|*2- 

2at/5 
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u*w; 

mi 

<h. 
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0 

& 

0 

0 

u\ 

Ui 

Ui 

m± 

Ujl 

0 

0 

0 

0 

U\ 

U 1 

C/1 

UfUi 

Ui 

0 

0 

0 

0 

£2. 

Vl 

Ui 

1/1 

(B.10) 

♦ 
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1 

mm:-- 

'  ■ 

-H2H1 

Hi 

Hi 

It 

in 

■■Mill 

u\ 

u  1 

U  i 

dp  u] 

Bp 

2  Hl+*P- 

dp 

dp 

dp 

dp 

U\ 

du2 

U 1  dU2 

dUi 

dU5 

dU6 

dUj 

u'dul 

Uzw5 

rr  dh0 

Ulw; 

U-JthsL 

*du7 

H1H1 

H I 

0 
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Hi 
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u\ 

U  i 
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Hl 
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0 

Hi 
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u\ 

U  i 

Ui 

H1H1 

Hi 
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0 

0 

Hi 

V\ 

Ui 

Ul 

(B.ll) 

Since  Wjj  is  not  required  for  either  stability  or  spatial  accuracy  this  matrix  is  not  computed. 
This  reduces  the  scheme's  time  accuracy  to  first  order  (34). 
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